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ABSTRACT 



This thesis presents the setup and results of numerical calculation of hadronic 
matrix elements of AS = 1 weak operators, with the aim of studying the AI =1/2 
rule and direct CP violation. Such study provides a useful comparison of the Stan- 
dard Model predictions with experiment. We work within the framework of (partially 
quenched) Lattice QCD with the staggered (Kogut-Susskind) version of fermions. We 
have achieved a reasonable statistical accuracy in calculating all matrix elements in 
question. Based on these results, we find that the AI = 1/2 channel in K — > nn 
decays is enhanced compared to the AI = 3/2 channel. This is consistent with 
experiment, although the exact amount of enhancement is subject to considerable 
systematic errors. I also discuss the difficulties encountered in attempting to calcu- 
late the direct CP violation parameter e' / e and our approximate solution leading to a 
crude estimate. In a related study we nonperturbatively compute bilinear renormal- 
ization factors Z$ and Zp, which are needed for our estimate of e'/e and also allow 
to calculate the light quark masses. For the strange quark mass in NDR MS scheme 
at 2 GeV in the quenched approximation we find: m s = 102.9 ±5.1 MeV. 
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a related study we nonperturbatively compute bilinear renormalization factors Z s 
and Zp, which are needed for our estimate of e'/e and also allow to calculate the 
light quark masses. For the strange quark mass in NDR MS scheme at 2 GeV in the 
quenched approximation we find: m s = 102.9 ±5.1 MeV. 
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CHAPTER 1 



INTRODUCTION 

1.1 Overview 

One of the main directions of efforts in both experimental and theoretical particle 
physics is obtaining reliable information about the fundamental free parameters of the 
Standard Model (SM). These parameters include the masses of quarks, leptons, W and 
Z gauge bosons, and the Higgs particle; coupling constants of electroweak and strong 
interactions; and the elements of the Cabibbo-Kobayashi-Maskawa (CKM) matrix. 
Without the exact knowledge of all of these parameters, the SM lacks the full predic- 
tive power. Such knowledge would also allow to perform various consistency checks 
based on results of independent experiments. In addition, there is another important 
reason for knowing the SM parameters: many alternative theories of fundamental 
particle interactions have been proposed. Among such theories are various extensions 
of the SM, as well as Supersymmetric Grand Unified Theories (SUSY GUTs), in 
which the SM itself is a low-energy limit of more fundamental and unifying theories. 
In order to check if extensions of SM are really necessary and to check if any of the 
proposed alternatives to the SM are indeed valid descriptions of physical reality, it is 
necessary to know the parameters of the SM. 
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In many cases obtaining the SM parameters from experiment is difficult or impos- 
sible due to the presence of non-perturbative effects. For example, in order to extract 
certain CKM matrix elements, very often it is necessary to calculate matrix elements 
(MEs) of weak interaction operators between strongly interacting states (hadrons). 
In these cases the weak and strong interactions are tightly intertwined; hence there is 
an obligatory non-perturbative component in computing MEs, since the energy scales 
involved fall within the nonperturbative region of Quantum Chromodynamics (QCD), 
the theory of strong interactions. In calculating this nonperturbative component, as 
in other aspects of nonperturbative calculations in particle physics, Lattice Gauge 
Theory (LGT) plays an increasingly significant role. It is a first-principles method, 
and although several types of errors are inevitably introduced by lattice computations, 
rapid advances in available computational power as well as algorithmic techniques are 
allowing for steadily better control of such errors. 

This thesis deals with several cases of application of Lattice QCD (LQCD) to 
calculating hadronic matrix elements of weak operators with the purpose of bridging 
the gap between theory and experiment and fixing the parameters of the SM. In 
addition, such calculations are aimed at testing nonperturbative QCD calculations 
against experiment so that the developed machinery of lattice calculations can be 
employed with confidence in other useful applications. 

This work addresses the rich phenomenology of K° — > tttt decays. One of the 
long-standing puzzles is the "A/ = 1/2 rule", which is the observation that the 
transition channel with isospin changing by 1/2 is enhanced 22 times with respect to 
transitions with isospin changing by 3/2. Weak interactions and the short-distance 
part of strong interactions are not enough to explain this phenomenon. In order 
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to study the long-distance part of strong interactions (namely, to calculate hadronic 
matrix elements of the operators appearing in the effective weak Hamiltonian) a non- 
perturbative method has to be employed. In the work of this thesis, LQCD is used 
to this purpose. The aim is to explain the AI = 1/2 rule and to check the validity 
of several SM assertions used in obtaining the explanation, as well as to test the 
technology of lattice calculations. 

In addition, I address the somewhat related issue of e', the direct CP violation 
parameter in the neutral kaon system. Experimentally, the value of e'/e is measured 
to be non-zero. Resolving a 6 years old controversy, the Fermilab KTeV group U 
has recently announced their new measurement of Re {e'/e) at (28.0 ± 4.1) ■ 10 -4 . 
This is consistent with the latest reports from CERN group NA48 |2j of Re [e 1 j e) = 
(18.5 ± 7.3) ■ 1(T 4 as well as with the value of Re {e'/e) = (23.0 ± 6.5) ■ 10" 4 measured 
by the CERN's previous experiment NA31 ||. Both Fermilab and CERN groups are 
currently analyzing their data, so enhanced statistics will soon allow to decrease the 
above errors. In addition, an experiment at Frascati is soon expected to measure the 
value of Re {e'/e) by a different technique. 

It is clearly very interesting to compare the theoretical prediction for e'/e with 
the experimental result. Such a comparison will directly show if the minimal version 
of the SM correctly accounts for direct CP violation. Estimating e'/e is currently 
an area of considerable effort in the theoretical physics community ||. According 
to some current estimates ||, it appears unlikely that the SM description of direct 
CP violation is correct, so that some extension or modification of the SM may be 
necessary. The weakest part of such estimates is the lack of knowledge of MEs of 
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certain weak operators [H], which require nonperturbative treatment. In this work, 
LQCD is used for calculation of the relevant MEs. 

The core of our work is calculation of matrix elements (7r7rlOjl.fr ) for all basis 
operators entering the effective weak Hamiltonian (introduced in Sec. |2|), within the 
framework of LQCD with staggered (Kogut-Susskind) fermions. Then the theoretical 
estimates for both AI =1/2 rule and e'/e can be made as mentioned above. 

Statistically significant results have been obtained for MEs of the basis opera- 
tors. Our results for AI = 1/2 rule (described in Chapter. |J) confirm a significant 
enhancement of AI = 1/2 channel compared to AI = 3/2 one. The exact ratio of 



isospin amplitudes (Sec. fO|) is subject to considerable systematic uncertainties, but 
is consistent with the experimental value. 

In principle, the MEs calculated in this work are sufficient to determine e'/e . 
However, currently the failure of lattice perturbation theory in operator matching in- 
troduces a large error in the end results. To reduce this error, nonperturbative opera- 
tor matching has to be performed. Such matching would allow to use the (statistically 
strong) results of this work to their fullest. To give a crude estimate, we have currently 



adopted a partial nonperturbative renormalization procedure (Sec. |5.3|) , combining 
perturbative mixing with nonperturbative determination of bilinear renormalization 
constants Z$ and Zp. This provides an admittedly rough way to do the operator 
matching. The estimates of e'/ e obtained in this manner come as a surprise: the sign 



of e'/e is negative (see Sec. |5.4| for numbers). This result contradicts the experimental 
findings. Of course, the current systematic errors preclude any definite conclusion, 
but if this situation persists in the future when such errors are reduced, the Minimal 
SM description of CP violation might need a revision. 
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Light quark masses are among the poorer known parameters of the SM. Yet, 
knowing their values is important, for example for constraining SUSY GUTs, where 
an independent prediction for these masses can be made. LQCD is one of the best- 
suited, direct techniques for extracting the light quark masses J7|. The major obstacle 
of lattice quark mass estimation has been the poor suitability of the perturbation 
theory for matching lattice and continuum operators. In this work such matching 
is done nonperturbatively, which provides better estimates of the light quark masses 
than previously available. For the strange quark mass in NDR MS scheme at 2 GeV, 
in quenched QCD we obtain: 

m s = 102.9 ±5.1 MeV . 
1.2 Comparison To Previous Work 

A variety of phenomeno logical methods such as 1/N C expansion [[§] and chiral 
quark models |J have been used for calculation of hadronic MEs relevant for both 
e'/e and AI = 1/2 rule. Without underestimating the strengths of such methods, 
it may be fair to say that in prospect LQCD is one of the best choices for such 
calculations (of course, assuming LQCD calculations are computationally feasible). 
This is due to the fact that LQCD calculations are intrinsically reliable: they are 
based on first principles in the sense that, by design, no unjustifiable approximations 
or assumptions are introduced. At the same time, it is also true that in many cases 
this strength of the LQCD method lies mainly in the future, since, as will become 
clear from this thesis, sometimes the uncertainties associated with current lattice 
calculations are overwhelmingly large. Most of these uncertainties will very likely be 
decreased in the future. 
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Concerning previous lattice calculations of the MEs discussed in this thesis, there 
have been several attempts, but they fell short of desired accuracy because of technical 
difficulties and/or insufficient statistics (due to limited computational power). Works 
reported in Refs. |I| |TT], [T2|, [13], |TJJ and reviews in Refs. [[TJ, [Tj| discuss attempts of 
several lattice groups to compute AI = 1/2 rule on the lattice, and explain various 
encountered difficulties. In addition, several groups [|17], [L2|, Q have studied MEs 
(n + 7i°\Oi\K + ) with reasonable precision; however, these MEs describe only AI = 3/2, 
but not AI = 1/2 transition. In the present simulation, the statistics is finally under 
control for both A J = 1/2 and A J = 3/2 amplitudes. 

A previous attempt |10| to compute MEs relevant for e' / e on the lattice with stag- 
gered fermions did not take into account operator matching. In this work we repeat 
this calculation with better statistics and better investigation of systematic uncer- 
tainties. In particular, we take into account operator matching by using a partially 
nonperturbative renormalization procedure (Sec. |5.3j ). As already mentioned, and as 
will be discussed in further detail in Chapter |5], operator matching is currently the 
source of the biggest error in computing e'/e. 

Recently, the RIKEN-BNL-Columbia group |18|] has reported the results of their 
first study of e'/e on the lattice using domain wall fermions. Compared to the stag- 
gered fermions, the domain wall fermion technique has the advantage of preserving 
the chiral symmetry away from continuum limit without introducing extra fermion 
flavors. The disadvantage of this method is the need to simulate a lattice with an extra 
(5th) dimension. The RIKEN-BNL-Columbia group is investigating whether simula- 
tions with domain wall fermions are practically feasible, and whether the advantages 
of domain wall fermion technique is worth the increased demand for computational 
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power. In any case, it is interesting to compare results obtained with these different 
methods. According to the latest work of the RIKEN-BNL-Columbia group [1I8J , the 
calculations with domain wall fermions also produce negative sign of e'/e . We note 
that our results have higher statistical accuracy and come closer to the continuum 
limit {(3 = 6.2) than those of the RIKEN-BNL-Columbia group. On the other hand, 
the RBC group uses fully nonperturbative matching, which is definitely an advantage 
compared to our partial scheme. The consistency of both groups' results at (3 — 6.0 
is very encouraging, and gives more credibility to our partially nonperturbative pro- 
cedure. 

Concerning calculation of light quark masses, only recently such calculations have 



started to use nonperturbative matching. Recently, JLQCD collaboration |19| has 
done a study of light quark masses with staggered fermions using nonperturbative 
matching. The work in this thesis uses similar ensembles, but a different method: we 
deduce Z$ and Zp from considering the inverse averaged propagator, while JLQCD 
group obtains these two factors from setting renormalization conditions on appropri- 
ate amputated Green's functions. There is also a small technical difference in the 
calculational setup. Despite these differences, our results are consistent: JLQCD ob- 
tains 106.0 ± 7.1 MeV, while we obtain 102.9 ± 5.1 MeV for the strange quark mass 
in NDR MS scheme at 2 GeV in quenched QCD. 



In addition, some work has recenetly been reported [20] on nonperturbative quark 
masses with Wilson fermions. This is a completely independent calculation since 
Wilson fermions are conceptually different from staggered fermions in many respects. 
It is very encouraging to see that the staggered and Wilson results are consistent. 
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Phenomenological methods, such as QCD sum rules, can also be applied for de- 
ducing quark masses [2T|. However, they possess several major disadvantages when 



compared to lattice calculations, as explained in Ref. [^2 
1.3 Organization Of This Thesis 

This thesis is structured as follows. In Chapter § we explain the framework 
in which we work: effective weak theory obtained by Operator Product Expansion 
(OPE) and Renormalization Group (RG) evolution. We also discuss the exact con- 
nection between the MEs that we strive to compute and the physical quantities of 
interest, as well as some other questions related to continuum physics. 

Chapter |3| delves into the theory and practice of lattice calculations. It starts by 



an overview (Sec. |3J]) of the generic setup of lattice calculations of hadron masses 
and MEs, including the most commonly introduced errors. This is meant to be not 
a systematic introduction, but rather a quick refresher. Staggered fermions are dis- 



cussed in a little more detail in Sec |3.1.5| , although the basic facts are only stated, not 
proven. The interested reader will readily find the missing information in available 
texts and review articles. Then, concentrating on lattice calculations relevant for 



this thesis, Sec. [O] discusses some matters of principle in such calculations. Sec. [O 
contains detailed explanation of the techniques and tricks we have used to implement 
the above strategies on the numerical level. Some of the techniques of dealing with 
staggered fermions are rather tedious, so I give explicit formulae for reference pur- 
poses. (Some other explicit expressions which may be useful to anyone attempting to 



perform similar calculations, can be found in the Appendix). Sec. [O] also contains 
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the parameters of our simulations. I close the chapter with a brief mention of some 
computational issues (Sec. |3.4| ). 

Turning to results of our calculations, Chapter [| presents the calculated AI =1/2 
and AI = 3/2 amplitudes and their ratio, with a discussion of the systematic errors. 
Chapter [5] is devoted to operator matching and e'/ e . Sec ^TT| explains how the op- 
erator matching problem together with other systematic errors precludes a reliable 



calculation of e'/e. Sec. |5.2| concerns calculation of Z$ and Zp bilinear renormaliza- 
tion factors nonperturbatively, including the detailed setup of lattice calculations for 
this case as well as lattice parameters used. Apart from results for Z$ and Zp, given in 



Sec 5.2.3, I present our results for the light quark masses obtained with nonperturba- 



tive matching in Sec 5.2.4. Sec 5.3 describes the partially nonperturbative procedure 



for matching four-fermion operators relevant for estimating e'/e , and Sec. |5.4| con- 
tains the results for e'/e and a discussion of the errors in its determination. A brief 
summary of results for e' / e is given in Sec. |5.5| . Chapter |5] contains conclusions. 
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CHAPTER 2 



THEORETICAL FRAMEWORK 

2.1 Effective Weak Theory 

The standard approach in applying the SM to topics mentioned in the introduction 
is to use the Operator Product Expansion (OPE) at the Mjy scale, that is to find a 
suitable basis of four-fermion operators which give effectively the same interactions as 
the charged- current weak interactions for a given problem of interest (see Ref. |§ and 
references therein). This is done by matching the strengths of interaction vertices 
of the basis operators with that obtained by a W exchange. The effective weak 
Hamiltonian for nonleptonic kaon decays is given as a linear combination of basis 
operators 0\ and O2, defined as follows: 

O x = (s Q 7 M (l - 7 5 )u ( g)(%7 M (l - 75)4,) 

2 = (Sa7/i(l - 75)M a )(%7 M (l - 75)40 

To compute physical amplitudes, one needs to compute matrix elements of this Hamil- 
tonian between the suitable hadronic states (in this work, kaon and two-pion states). 

Note that the most natural energy scale \x for lattice calculations is around or less 
than 1/a, which is 2 - 4 GeV in current simulations. This is much less than the scale 
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Mw of the above effective Hamiltonian. It follows that if the QCD corrections are 
to be included in the effective theory description, their contribution from one-loop 
diagrams (which are 0(a s log(m^//i 2 ))) can be sizable. Higher-order corrections are 
0([a s log^rri^Y / fi 2 )] n ) , which should not be neglected. A well-known method of re- 
summation of these "leading logarithm" terms involves using Renormalization Group 
(RG) equations to evolve the theory down the energy scale to a region more accessible 
for current lattice calculations (2-4 GeV). Then the resulting theory is valid up to 
terms 0(a s [a s log(m^//i 2 )] n ), and the method is called RG Improved Perturbative 
Expansion. To achieve an even higher accuracy, as is quite necessary in the case of 
kaon nonleptonic decays, next-to-leading order logarithms can be resummed by using 
one-loop matching at the Mw threshold and using RG equations with anomalous 
dimension matrix computed to two loops. 

In the process of RG evolution heavy quarks are integrated out by matching the 
theories with and without the given quark flavour at a threshold of the order of the 
given quark's mass. For detailed discussion of this procedure refer to Ref. 0. Due 
to elimination of heavy quarks from the effective theory and due to operator mixing 
the basis of operators is enlarged to 10. The effective Hamiltonian takes the following 
form: 



where Gp is the Fermi constant, z% and are Wilson coefficients (at two- loop order), 
r = —V t dV t */V u dV* s , V is the CKM mixing matrix, and Oi is a basis of four-fermions 
operators defined as follows: 
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(2.1) 



i=l 



o 1 



(Social ~ l5)up)(upY(l - J 5 )d a ) 



(2.2) 



2 



(Sa7/i(l - 75)«a)(W/37 M ( 1 - 75)0?/3) 



(2.3) 
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(s a ^(l - 


75K)^(g/37 M (i + 75)^) 


(2.6) 


o 6 = 


(s a ^(l - 


75)^) E(^7 M (i + 7 5 )ga) 
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(2.7) 


7 = 




- 75K) J! e q (q/rf(l + 75)^/3) 
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(2.8) 
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- 75)^/3) e 9(9/37^(l + 75)g Q ) 
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(2.9) 


o 9 = 


2(^(1 


- 75)^) I] e 9 (^7 /1 (l - 75)^) 
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(2.10) 


o w = 


2(^(1 


9 


(2.11) 



Here a and /3 are color indices (summation over which is implied), e q are quark electric 
charges, and summation over q is done over all light quarks. This description is valid 
up to 0(a 2 s [a s \og{m% / / n 2 )] n ) terms. 

2.1.1 Treatment Of Charm Quark 

A natural question at this point is whether one should consider the charm quark 
as "light" or "heavy" . In the first case, it is to be included in the sum over q in the 
above equations, and some additional operators containing the charm quark should 
be considered; while in the case it is considered "heavy" it is integrated out, and only 
u, d and s quarks are left in the theory. Since the scales we consider are of the same 
order as m c , the choice is ambiguous. The physical results are not expected to depend 
on this choice, disregarding the 0(a^[a s log(m^//x 2 )] n ) terms and higher-dimensional 
operators (proportional to l/m c ). 

In practice, when performing lattice calculations it is more convenient to work 
in the three-quark effective theory. In this framework the set of operators is smaller 
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than it would be in the four-quark theory, and the problems of discretizing the charm 
quark field are not present. Such approach is adopted in the present work. 

2.2 Definitions 



Here we define some quantities of interest. The AI =1/2 and AJ = 3/2 transition 
amplitudes are defined as follows, according to the isospin of the final two-pion state: 

A 0>2 = ((tttt) i=0 , 2 \H w \K°) e"^ 2 , (2.12) 

where 5o, 2 are the final state interaction phases of the two channels. Experimentally 

to = ReA /ReA 2 ~ 22 . (2.13) 

One goal of our work is to check if the theory gives the same ratio. Also, we would 
like to verify if predictions for ReA and ReA 2 separately agree with experiment. 

Another direction of our work is calculating the size of direct CP violation. Direct 
CP violation refers to the fact that the neutral kaon state K 2 , defined as an odd 
CP eigenstate, has a non-zero (though tiny) probability to decay into two-pion state, 
which is CP-even. This is to be contrasted with the indirect CP violation, a bigger 
and better-studied phenomenon arising from the fact that the eigenstates of the weak 
Hamiltonian (K s and K L ) are mixtures of CP eigenstates (Ki and K 2 ). Indirect CP 
violation is characterized by parameter e, which is defined as 

_ A(K L -> tttt) 
£ ~ A(K S -> tttt) ^- Mj 

and experimentally measured at 2.26 • 1CT 3 . Direct CP violation is characterized by 
parameter e', defined as 

A(Ki -> tttt) 
13 



As mentioned in the introduction, two latest experiments report the measured ratio 
Re {e'/e ) at (28.0 ± 4.1) • 1(T 4 and (18.5 ± 7.3) • 1(T 4 . 

It can be shown |23|, |24| that in the Standard Model e' can be computed in terms 



of imaginary parts of the isospin amplitudes defined above: 

, Iim4n — uj\m.Ao s( /,,. x \ , 

e = ^ e J ^/ 2+52 - 5 °). 2.16 

y/2 u ReA 

Expressing the amplitudes in terms of the basis operators of the effective Hamiltonian, 
we find that the experimentally measured ratio of direct to indirect CP violation is 
given by 

Re(^)= ^ ImAt [n -^n 2 ], (2.17) 

where 

n = £yi((7r7r) J=0 |of \K°)(i-n v+v ,) (2.18) 

i 

n 2 = Y,yi(t™)i=*\ T ] \ KQ ) ( 2 - 19 ) 

with lm\ t = lmVtdV t *, and where f^+j/ ~ 0.25 ± 0.05 takes into account the effect 
of isospin breaking (since m u ^ md). of and O t - are isospin and 2 parts of the 
basis operators. Their expressions are given in the Appendix for completeness. The 
matrix elements in the last two equations contain long-distance dynamics of strong 
interactions, and therefore a non-perturbative calculation, such as the one considered 
in this thesis, is necessary to compute e'/e . 
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CHAPTER 3 



LATTICE SIMULATIONS 



As explained in the preceding chapter, crucial knowledge about long-distance dy- 
namics is contained in matrix elements (7nr\Oi\K°) of weak operators from the basis 



in Eqs. (|2.2| )- (]2TTI ). In this chapter we review the basics of lattice calculations, dis- 
cuss the setup and techniques of calculation of matrix elements in question, and give 
detailed information about parameters used in the simulation. 

3.1 Lattice Background 

Let me start with a quick overview of general facts about Lattice Gauge Theory 
(LGT) calculations. More details are available to the interested reader in numerous 



textbooks and review articles |25, ^6|, ^8], |29[. LGT is a systematic, first-principles 



method for numerical solution of various theories which are hard or impossible to 
solve analytically. In particular, one of the most often mentioned strengths of the 
LGT technique is its ability to deal with non-perturbative aspects of theories such as 
QCD. 

LGT has been successfully applied and has made a number of important contribu- 
tions in various areas such as calculation of matrix elements, quark masses, glueballs, 
heavy flavour physics, finite temperature theories to mention a few. The impact of 
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LGT on particle physics phenomenology is likely to expand even further in the fu- 
ture, since most of the uncertainties of LGT simulations can be improved without 
limit with increasing computational resources. In this work I am going to concentrate 
on application of LGT to calculating matrix elements of weak operators as well as 
light quark masses in QCD. 

Computer simulations are made possible by discretizing Eucledian space-time on 
a four-dimensional hypercubic lattice of finite extent in all directions. This certainly 
introduces various errors compared to the continuum description, for example, viola- 
tion of Eucledian invariance and introduction of ultraviolet and infrared momentum 
cutoffs. In practice, the discretized version of the theory is formulated in such a 
way as to make possible the extrapolation to continuum limit by repeating identi- 
cal calculations at progressively smaller lattice spacing values (denoted a). Since 
one wants to keep the physical size of the "box" fixed, the continuum extrapolation 
means increasing the number of lattice sites, and therefore making the computations 
progressively more computer-intensive as one approaches the limit of a = 0. One also 
has to check that the physical size of the box is large enough to contain the system in 
question. This can be established by checking that results do not change in repeated 
calculations at bigger volumes. Most often periodic boundary conditions are used. 

3.1.1 Lattice QCD Lagrangian 

Lattice QCD simulations are done in the framework of Eucledian functional in- 
tegral formulation of the theory. As is well-known, the continuum QCD Lagrangian 
equals 
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The partition function can be written as 

Z = JVU Vi)V$ e -s a [U]-s q [^m = J vU~[[det(M q [U})e- s ° [u] , (3.2) 

where det(M q [U]) = det($> + m q ). 

In the original formulation of Lattice QCD by K. Wilson |30| the quark fields are 
defined on each lattice site, while the gluon fields "live" on links connecting these 
sites. Namely, consider a gauge link, an SU(3) matrix variable corresponding to the 
path-ordered line integral along one link in the continuum: 

U„(x) = P exp(ig [ X dx' A c ^x')T c ) , (3.3) 

where T c are color group generators, x is a lattice site, g is the QCD coupling constant. 
In the lowest-order in a, these variables are associated with gauge fields as follows: 

U^x) = exp(-iagAfa + fia/2)T c ) . (3.4) 

The gauge links provide a convenient way to introduce gauge-invariance. Basically, 
all gauge-invariant quantities are either products of gauge links along a closed path in 
space-time, or expressions of the type V(x)f/(a;, 2/)?/^(2/), where U(x,y) is the product 
of gauge links along a path from y to x. 

The construction of lattice QCD Lagrangian is a standard textbook material. The 
gauge part of action equals 

3 g [U] = /?E EC 1 - ^ Re TrP M ,(x)) , (3.5) 

X fJ,<V " 

where (3 = 6/g 2 and P^ v (x) is so-called plaquette, defined as the product of 4 gauge 
links around the minimal closed path: 

P^(x) = Uy^U^x + ua)Ul(x + va + /ta)f/^(x + fia) . (3.6) 
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The fermionic part of action is 

S 9 f0,V>,£/]=£W> + m)V> , (3.7) 

X 

where p is a suitably discretized, Euclidean covariant derivative (to be defined later). 
Quantities of interest in lattice simulations can be extracted from correlation functions 
obtained by weighted average as follows: 

(C) = \ [ VU]Jdet (M q [U]) e~ s ° [u] C[U] . (3.8) 
Z J q 

Taking such integrals is an enormous task, since it involves integrating over 4L 4 SU (3) 
variables. In practice, Monte Carlo with importance sampling is used. A number of 
efficient algorithms have been developed For this purpose, the simplest of them being 
the Metropolis (heatbath) algorithm. We only mention the names of the algorithms 



used for this work in Section [3.3| however, the discussion of them can be found in 
many texts on lattice gauge theory and is beyond the scope of this work. The bottom 
line is that any algorithm produces a set ( "ensemble" ) of iV gauge configurations [/, 
at given (3, which are randomly sampled from probability distribution 

P oc ]Jdet(M q [U])e- S9[u] . (3.9) 

Correlation functions are then obtained by averaging over the ensemble: 

1 N 

(c) = Tr T,cm • (3-10) 

iv 1=1 

3.1.2 Quenched Approximation 

It should be mentioned that in practice calculation of fermionic determinants in 
Eq. ( |3.9|) is computationally very demanding. Only several years ago realistic lattice 
calculations have started to include it in simulations. The alternative is to replace it by 
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a constant. This is the so-called quenched (or valence) approximation. It corresponds 
to a continuum theory just like QCD but without any closed fermion loops. In 
other words, it is assumed that the effects of the sea quarks can be absorbed into 
renormalized charges. The error introduced by such an approximation is not known 
exactly, but in most cases Quenched QCD (QQCD) is expected to produce results 
within 10-15% of those obtained in QCD. This has been confirmed numerically for a 
number of quantities. Thus, QQCD can be thought of as a convenient test platform 
for developing algorithms for many calculations, being very similar to QCD. Moreover, 
the results obtained in QQCD are sometimes valuable by themselves, especially in 
cases where it is interesting to know the answers even at the order of magnitude. 

In this work we have mostly used quenched approximation. To see if quenching in- 
troduces any serious error in the quantities of interest, we have also used one partially 
unquenched ensemble with 2 flavors of sea quarks of mass m = 0.01 at f3 — 5.7. 

3.1.3 Extraction Of Hadron Masses And Matrix Elements 

First consider correlation function 

(C(t)) = <0|T[Ot WO( o)]|0> , (3.11) 

where O is an operator that has the suitable quantum numbers to create a hadron. 
For example, in order to create n + at rest, the operator 0(t) = J2 x u(t,x)-f 5 d(t,x) 
can be used. After insertion of complete set of states and assuming t > we obtain: 

(C(t)) = (0|Ot( ) e -«O(0)|0> = £ |<n|O|0)| ^ . (3.12) 

If we consider operators of zero spatial momentum, then the lightest particle for which 
(n|O|0) is non-vanishing dominates the latter expression for large t. In cases of our 
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interest it is always the pion or kaon []. Making an exponential fit to the correlator 
C(t) data, it is possible to extract the mass of the hadron. In addition, it is possible 
to extract the matrix element (7r + |w75<i|0). If we consider operator u^^^d, we can 
extract the pseudoscalar decay constant f n (see Sec. |4.1.2|) . 

In order to show how the calculation of C(t) is done, let us rewrite Eq. ( |3.11| ) as 
follows: 

(C(t)) = $>|T[3 Qi (x,t) r y u aj (x,t)u pk (y,0) T kl ^(y,0)|0) , (3.13) 

where the summation is done both over color and spin indices and T are any matrices 
specifying the spin structure of the operators. In terms of the propagator G(m;x,y) 
of a quark of mass m defined by 

G?f(m;x,y) = (0\T[^(x)^{y)]\0) (3.14) 

the expression for (C(t)) can be rewritten as follows^: 

(C(t)) = E(^(x,t;y,0)r fe/ Gf(y,0;x,t)r,) 

= ^(Tr[G(x, t- y, 0) V G(y, 0; x, t) T}) , (3.15) 

where the trace is taken over both color and spin indices and the brackets imply 
averaging over gauge configurations. 

For our purposes it will be necessary to calculate three-point functions for matrix 
elements of the type {n\0 k-k\K) , and here I would like to mention the setup of such 

calculations. Suppose the kaon source is at to = 0, the operators are inserted at an 

-'^Note that in our simulations we have always taken m s = m c i, so that the strange and down 
quarks are indistinguishable in terms of computations, and therefore kaons are indistinguishable 
from pions. 
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The corresponding expressions for staggered fermions are more complicated, and they are given 



below in Sec. 3.3, along with detailed explanation of the meson sources and sinks that we have used. 
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t 



T 



Figure 3.1: The general setup of the simulation. An "eight" contraction is shown 
for convenience. The kaon source is at the timeslice 0, while the pion sink is at the 
timeslice T. The operator is inserted at an adjustable time t. The result of this 
contraction is proportional to the product of two exponentials shown in the figure. 



adjustable time £, and the pion sink is located at the time T (see Fig. p.l|) , where T 
is sufficiently large. Consider the correlation function 

(C(t,T)) = (0\O v (T)O Kv (t)O K (0)\0) . (3.16) 

Inserting two complete sets of states as before, we obtain 

p—E ni (T—t) p—E„ 2 t 

(C(t,T))= ^(OjO.COK) {ni \o Kw (t)\n 2 )——-(n 2 \O K (0)\0). (3.17) 

As before, the lightest states (in this case, kaon and pion) dominate the expression 
at large t, so we obtain: 

(C(t, T)> - Ze- m * T (n\0 K7T \K) , (3.18) 
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where yfZ = {^O^tt) 2r ^ y i and we have assumed that both kaon and pion are created 
at rest and have equal masses. Note that the last expression does not depend on t, so 
we expect to see a plateau on the plot of (C(t, T)) vs. t. This plateau is our working 
region (see Fig. |3.4| ). It is present only if T is large enough. In order to know exactly 
how large T should be in order to be free from excited states contamination and 
wraparound effects, we have done some checks (see discussion in Sec. |3.3.4|) . Explicit 
expressions for (C(t,T)) in terms of propagators are given below in Section |3.3| . 
In addition, we will need to compute the ratio of matrix elements 

mm ( 3 19) 

in which the factors \fZe~ mt cancel and one expects to see a plateau for large enough t. 

Note that all quantities extracted from the lattice are in lattice units. In order to 
convert to physical units, each quantity should be multiplied by a~ p , where p is the 
mass dimension of the quantity. The lattice spacing uniquely corresponds to (3. In 
particular, in the perturbative scaling region, when the lattice spacing is small enough, 
the following relationship exists: 



16** ; p U/w 



a{(5) = a exp j (3.20) 



'MS' 

where (3q = 11 — §-/V/, f3\ = 102 — yiVy and Nf is the number of sea quark flavors. 
The coupling constant can be obtained from the average plaquette value at a 
given p (see Eqs. ( |3.62| ) and ( |3.63| )). The overall scale ao can be obtained by dividing 
one chosen quantity in lattice units by the experimental value in physical units. For 
example, the mass of the p meson can be calculated on the lattice, and the lattice 
spacing can be set by 

(rapO)latt 



m p (GeV) 
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Exactly how the lattice scale is determined varies from computation to computation. 
In order to take the continuum limit, calculations should be repeated with several 
values of (3 in the region of small enough a. 

3.1.4 Discretizing Fermions 

So far we have not discussed the fermion part of the Lattice QCD action, in 
particular the discretized version of the operator Jp. One straightforward (naive) 
scheme is to write the fermion action as 

$N = o ^(xhtilU^x^ix + ap) - Ul(x - aft)ip(x - afi)} + ^ mtp(x)^p(x) , (3.21) 

^ X,fJ, X 

where 7 M are the Eucledian Dirac matrices, which are hermitian and satisfy {7^, 7^} = 
25^. This form of action satisfies the important gauge invariance requirement and 
seems to correspond to the QCD action in continuum limit. However, a well-known 
problem with the above naive action lies in the following. According to this action, 
the free fermion propagator is 

G(p) = . (3.22) 

1 J2f, sm(p / ,a)7 At + m 

This propagator implies 16 poles at p e {0,7r/a} 4 (in the limit of vanishing fermion 
mass m), which includes not only the physical pole in the vicinity of 0, but also 15 
spurious poles corresponding to other species of fermions (so-called doublers). Another 
important feature of the naive action is that the doublers come in pairs of opposite 
chirality, so that the net chirality is always zero. In particular, the chiral anomaly is 
completely absent |3l| . Using the naive action, it is impossible to define any theory 



with chiral fermions, such as the electroweak theory. Moreover, the naive action is 
unsuitable for describing QCD, since the chiral SU(3) x SU(3) symmetry, essential 
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to the theory, is reduced to vector SU(3), and important axial Ward identities are 
missing. 

Of course, the naive fermion action above is not the only possible action with 
needed continuum limit. Much work has been done in the direction of devising a 
fermion discretization scheme that would be free from doublers and, more importantly, 
have the correct chiral properties. In such work, one has to be aware of the Nielsen- 
Ninomiya theorem [p2], stating that the fermion doubling problem is unavoidable 
for any local, anti-hermitian discretization with periodic boundary conditions that 
preserves translational invariance. 

Currently most popular approaches which in some way deal with the "doubling 
problem" include Wilson fermions, domain wall fermions and staggered fermions. In 



the approach of Wilson fermions [g^] an extra term is added to the action which 
changes the propagator in such a way as to make the effective physical mass of the 
15 extra doublers infinite in the continuum limit; therefore they effectively decouple. 
Unfortunately this extra term completely breaks the chiral symmetry, though it may 
be recovered in the continuum limit. The approach of domain wall fermions is to 
overcome the Ninomiya-Nielsen theorem by introducing the 5th dimension of the 
lattice. This approach preserves chiral symmetry without extra doublers. However, 
the computational cost of dealing with the 5th dimension is high. The Wilson fermions 
and domain wall approaches have been described elsewhere in detail. In the following 



I concentrate on the approach of staggered (Kogut-Susskind) fermions fl34 |, since it 
was assumed throughout this work. 
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3.1.5 Basic Facts About Staggered Fermions 

In the staggered fermion approach, the extra doublers are kept. However, their 
number is reduced from 16 to 4 (by using the spin-diagonalization procedure, see be- 
low), and in such a way as to obtain nice chiral properties even away from continuum. 

The spin-diagonalization procedure is done by changing variables from x to 



4>(x) = T x x(x) , 
where 4x4 matrices T x satisfy the following property: 

with rjulx) being diagonal unitary matrices. One possible realization is to use 

r« = it'iTifiT ■ (3.23) 

With this choice the phases rj^ are equal to 

7fc(x) = (_1)^+-+-m-i . (3.24) 
It is easy to show that the naive action in terms of new variables is 
S = ^^^(x)[x(x)f/ M (x)x(x + a/i)-x(x + a/i)^(x)x(x)]+^ m x(a ; )x(a ; ) • (3-25) 

X,fJ, X 

Since this action does not mix four components of x, they contain redundant informa- 
tion, and therefore three of them may be dropped. From now on we will understand 
that in Eq. ( p.25|) x( x ) an d xi x ) are single- component Grassman variable fields. Note 



that in the kinetic term the fields at odd sites are coupled only to fields at even sites. 
Therefore, in addition to the symmetries of charge conjugation, parity, shifts and dis- 
crete rotations, the staggered action preserves the U 0< a(l) x U even (l) symmetry in the 
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massless limit, which corresponds to a combination of fermion number conservation 
and U(l) chiral symmetry in the continuum. 

The spin-diagonalization procedure allows one to decrease the number of fermion 
doublers by 4 times. The remaining 4 species can be identified as independent 4 flavors 
of physical fermions in the continuum as follows. In order to "make contact" with 
physical Dirac spinor fields, the space-time is divided into hypercubes, each spanning 
2 sites in each of the 4 dimensions. The field x at the 16 sites of each hypercube is 
interpreted as various spinor components of the "physical" quark field. In this way, 
one introduces quark fields that "live" on a coarse lattice with spacing 2a and have 16 
components, which correspond to 4 Dirac (spin) indices and 4 indices corresponding 
to independent, identical flavors of fermions in the continuum limit. Mathematically, 
one defines 

f a (y) = W T A;aaX(y,A) , (3.26) 
8 A 

T a {y) = \j2x(y,A)rl taa . (3.27) 

A 

In this equation \ and X are labeled by hypercube coordinate y and coordinate inside 
the hypercube A G {0, l} 4 , so that the total coordinate of a given hypercube corner is 
x = (2y+Aa). The sum runs over the 16 corners of a hypercube. Spin (flavour) indices 
are labeled by Greek (Roman) letters. By using the properties iTr[r^r B ] = 5ab and 
\ Ha ^A-,b/3^A;aa = 5 a p5 a b-, it is possible to show that the free action takes the form 

S = 16^m(5(y)l®l g (y)) + 16^^g(y)( 7/1 «)l)A /1? (y) 
y y m 

+ 16a£E9fo)(75®U5)W2/)> (3.28) 
y m 

where ^ = 7J and the direct product 7 <8> £ means that the 7 matrix acts on Dirac 
(spin) indices and £ acts on the flavour indices. and 5^ are the lattice versions of 
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first and second derivative: 



\q(y) = ^{q{y + A) - g(y - A)) — > d^( y ) , (3.29) 



6 Av) = ^(9(y + A) + 9(y - A) - 2g(y)) — 0jg(j/) . (3.30) 



As already mentioned, this action defines 4 flavors of fermions which are intermixed 
by the last term of the action. However, this term is proportional to a, and so in 
continuum limit the flavors completely decouple and become independent. 

To summarize, by introduction of quark fields on a coarse lattice, it is possible to 
find a correspondence between the single- component fermion field x an d spinor fields 
corresponding to 4 independent flavors in the continuum. The doubling problem of 
the naive fermion action is overcome, since now the propagator in momentum space 
has the form 

G( , ^ a I sin(2op^) fa ® 1) + a S M (cos(2ap M ) - 1)(7 5 g ^js) - ma{l g 1) 

X)„ sin 2 (ap M ) + m 2 a 2 

(3.31) 

and has only one pole (at p = 0) within the allowed range of momenta —ir /2a <p^< 
7r/2a, corresponding to the coarse lattice. The axial anomaly is non- vanishing and 
has the correct continuum limit. 

If it was not for the last term in Eq. ( |3.28|) , the action would possess a t/(4) vect0 r x 



^(4) ax iai symmetry in the massless limit. The last term breaks this symmetry down 
to U(l) x U(l), with the generators of symmetry transformations being V — l)q 
and A = g(7 5 <g> £5)^. If the axial symmetry is spontaneously broken, the associated 
Goldstone boson is flavour non-singlet, with the spin-flavour structure 75 <S> £5. This 
will be important in our calculations, wherein we evaluate matrix elements between 
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states of Goldstone bosons associated with spontaneous breaking of the chiral sym- 
metry. Thus, staggered fermions retain some basic chiral properties of the continuum 
theory. In particular, the mass is multiplicatively renormalized. A number of useful 
Ward identities have been derived in Ref. PB1 based on this property, which greatly 



simplifies our calculations by predicting the chiral behaviour of matrix elements. The 
chiral properties of staggered fermions are a significant advantage compared to the 
case of Wilson fermions, where one has to satisfy the Ward identities by fine-tuning 
parameters and considering an enlarged set of operators that can mix with the basic 
ones. This advantage arguably outweighs the disadvantage of dealing with flavour in 
addition to spin structure of operators, as well as intricate mixture of internal and 
space-time indices. 

In numerical simulations we use the action in Eq. ( |3.25| ). However, the interpre- 
tation of single- component fields \ in terms of spin-flavour is necessary if one wants 
to make contact with continuum physics. In practice, this interpretation is important 
when one constructs quark bilinears (or four-fermion operators) with given quantum 
numbers S (for spin) and F (for flavour). Then it is necessary to be able to write down 
such bilinear contractions in terms of the single-component field. For such purposes, 
we use the operators 



°(y) = E X(V, S)( 75 ® t F ) BA x{y, A) , (3.32) 

10 A,B 



where (75 (g> £f)ba stands for the matrices 75 ® £p in the following representation: 

(is®Zf) ba = ^Tr[r^r 5 r A rJ,] . (3.33) 
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Note that the expression in Eq. ( |3.32| ) is non-vanishing only for A and B such that 
+ B^ + Sfj, + = (0 mod 2). So one can rewrite it as 



1 



°(V) = T7:J2x(y,B)( 7s 0^ F ) BA x(y,A)6 B , A+A , (3.34) 

iD A,B 

where A M =2 + F^. Note that all operations involving A,B, S and F are under- 
stood modulo 2. Four-fermion operators can be constructed by multiplication of two 



bilinears. For explicit formulae used in calculations, see Sec. |3]3|. Note that in gen- 
eral, the flavour structure of any physically relevant operator does not matter since 
in the continuum all flavors are equivalent. However, dealing with Goldstone bosons 
is special since only the axial symmetry with quantum numbers 75 (g) £5 is exact in 
chiral limit, and this dictates a choice for the flavor structure of operators in many 
cases. The rules for the flavour structure in such cases were laid down in Ref. 



and are summarized below in Sec. 3.3 



So far, we have considered only the free fermion action. The description of 
fermions interacting with gauge degrees of freedom can be obtained by inserting 
a factor U(B, A) representing a product of gauge links along any path connecting 
the two sites. This makes each such expression gauge- invariant. The same is true 
for staggered fermion operators. (It is often advantageous, as in this work, to make 
all operators gauge-invariant. This restricts the number of operators allowed to mix 
with the original ones.) There is an ambiguity as to which path should be chosen. All 
choices are equivalent since they are supposed to give the same results in continuum 
limit. For technical purposes, we have chosen to use the product of gauge links along 
the shortest path, or, if this path is not unique, average over all the shortest paths. 



See more discussion on this in Sec. 3.3 
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3.1.6 Summary: The Skeleton Of A Lattice Calculation 

I have overviewed the basics of generic LGT calculations. Many important details 
regarding the calculations presented in this thesis can be found below in Sec. [373 



Here I give a brief outline of a lattice calculation from a procedural point of view. 

1. Discretize space-time on a 4D lattice of spacing a and a finite box size L. 

2. Work in the framework of Euclidean functional integration, with a suitably 
discretized QCD Lagrangian. 

3. Generate an ensemble of gauge configurations using Monte Carlo importance 
sampling. Each ensemble is characterized by a common (3 = Q/g 2 . Often, it is 
necessary to quench fermion loops in order to save computer time (quenched 
approximation) . 

4. Compute quark propagators from suitable sources, on each gauge configuration 
of an ensemble. 

5. Construct correlators of quark operators with quantum numbers appropriate 
for given hadrons. Extract physical observables (such as hadron masses and/or 
matrix elements) in lattice units. 

6. Average quantities over gauge configurations within a given ensemble. 

7. Set the lattice scale a~ l by using a suitable observable (for example, the mass 
of the p meson). Convert results to physical units by dividing by powers of a. 

8. Sometimes (as in this work) it is necessary to perform chiral extrapolation or 
interpolation: repeat calculations at several quark masses, then extrapolate or 
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interpolate in quark mass to obtain results at quark masses corresponding to 
physical hadron masses. Chiral perturbation theory is useful in such extrapo- 
lations for predicting the chiral behaviour of quantities of interest. 

9. Lattice operators should be matched to continuum ones, as discussed in Chap- 
ter 5. 

10. Repeat steps 3 — 9 for several values of (3 (i.e. different values of a), keeping 
the physical size of the box constant. Attempt to take continuum limit by 
extrapolating the quantities of interest to a = 0. 

3.1.7 Sources Of Error In Lattice Simulations 

Although Lattice QCD is a systematic, first-principles approach, certain statistical 
and systematic errors are present in current simulations. The hope of the future 
simulations is to decrease these errors. Here I mention the most common of them. 

1. The statistical error comes from a finite number of gauge configurations (finite 
sample size). It is present in virtually all lattice calculations. It can be decreased 
by increasing the number of configurations (the error behaves as 1/y/N). To 
estimate this error, a number of intricate error analysis techniques can be em- 
ployed, taking into account correlations between various quantities. 

2. The error due to finite cutoff effects exists in current simulations because the ex- 
trapolation to continuum limit is not reliable in many cases. This has to do with 
the fact that the calculations become increasingly computationally demanding 
as one approaches the a = limit. Apart from just increasing the computational 
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power, several approaches exist in reducing the finite cutoff effects, including 
improved and perfect actions. 

3. As mentioned, finite box size can introduce an error. Thus it should be checked 
that the box is large enough. 

4. Often, chiral errors are present. Since the inversion of the fermion matrix 
becomes computationally more expensive due to the phenomenon of "critical 
slowing down", one is forced to do the calculations at relatively high quark 
masses and extrapolate to lower ones. Such extrapolations introduce an error. 
In addition, sometimes (as in this work) the quark masses in a given hadron are 
taken equal, whereas they are not equal in the real world. 

5. Quenched approximation introduces an unknown error since it completely re- 
moves the effects of sea quarks. Theoretical as well as numerical estimates 
indicate that this error is not more than 10-15% for most quantities. However, 
one should be careful since this error can vary from one quantity to another. 
Ideally, one would like to perform full QCD calculations (including the effects 
of the sea quarks), but such calculations are presently 100 — 1000 times more 
expensive than quenched ones. 

6. There is an error associated with matching of lattice and continuum operators. 
I discuss it in detail in Chapter 5. 

This ends the overview of errors commonly encountered in lattice calculations. 
Discussion of these errors regarding the calculations presented in this thesis can be 
found in Chapters 4 and 5. 
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3.2 Theoretical Issues 



Before turning to detailed discussion of numerical simulation techniques, let us 
consider several issues of principle for our computations. 

3.2.1 Calculating (ttttIO^K ). 



As was shown by Martinelli and Testa [38], two-particle hadronic states are very 
difficult to directly construct on the lattice (as in any scheme using Euclidean space- 
time). We have to use an alternative procedure to calculate the matrix elements 
involving two pions in the out state. We use the method due to Bernard et al. 
which relies on chiral perturbation theory to relate (im\Oi\K ) to matrix elements 
involving one-particle states. That is, we use the following relationships which are 
obtained in the lowest-order of chiral symmetry breaking in QCD: 



(tt+tt-IO^ ) = ^~ m - 7 (3-35) 
(7r + \O i \K + ) = ip n .p K ) 7 -!^±I^d (3.36) 
(0\Oi\K°) = (m s -m d )5, (3.37) 



where / is the lowest-order pseudoscalar decay constant. The masses in the first 
of these formulae are the physical meson masses, while the quark masses and the 
momenta in the second and third formulae are meant to be from actual simulations 
on the lattice (performed with unphysical masses). 

The above relationships ignore higher order terms in the chiral expansion, most 
importantly the interaction between the two pions in the final state. Therefore this 
method suffers from a significant uncertainty. In principle, higher-order terms may 
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be computed. For example, Golterman and Leung [ JD| have computed one-loop cor 



rection for AI = 3/2 matrix elements in a closed form. However, the values of 
several parameters (so-called contact terms and the momentum cut-off) are currently 
unknown and therefore subject to guessing. Varying the unknown parameters in a 
reasonable range, one obtains that the first-order corrections are about 30% to 60%. 

In our work we have chosen to use strictly lowest-order relationships, remembering 
that this method is subject to considerable systematic error. (It is known that the 
AI = 1/2 amplitude is underestimated while the AI = 3/2 amplitude is overesti- 
mated by this procedure). If certain progress is made in the future to improve the 
quantitative estimates of higher-order chiral corrections, these developments can be 
easily incorporated into the present framework to produce more reliable answers. In 
other words, it will not be necessary to repeat the extensive lattice calculations in 
order to reevaluate the main results obtained in this work. 

3.2.2 Mixing With Lower-dimensional Operators. 



Eqs. ( |3.35| - j337D show that in order to recover {7nr\Oi\K°) , it is necessary to sub- 



tract from (ir + \Oi\K + ) the term proportional to 5, which corresponds to unphysical 
s <-» d mixing. This term is proportional to (0\Oi\K°). There are several ways to 
implement this subtraction in lattice simulations. We choose the method suggested 
in Refs. |B7), which is to subtract the operator 



O sub = (m d + m s )sd + (m d - m s )sj 5 d . (3.38) 

This operator, having dimension 4, has to be subtracted non-perturbatively from the 
basis operators (all having dimension 6). (A perturbative scheme would inevitably 
ignore some higher-order terms, which may be proportional to powers of a -1 and 
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therefore diverge in continuum limit.) Therefore, the mixing coefficients are deter- 
mined non-perturbatively by requiring that the subtracted (0\Oi\K°) vanish. Thus 
we have: 

2 2 

(n+n-lO^K ) = (n + \0, t - a i O sub \K + ) ■ ^ , (3.39) 

where an are found from 

= (0\O t -a t O sub \K°). (3.40) 



Kilcup and Sharpe |36[ derived a number of Ward identities which show that the lat- 



tice formulation of QCD with staggered fermions retains the essential chiral properties 



of the continuum theory. In consequence, as mentioned in Ref. [p7l, the procedure 



expressed in Eqs. ( p. 39 , 3. 40 ) is a faithful lattice implementation of the lowest-order 



chiral perturbation theory prescription, which reproduces Eqs. ( |3.35| - p37D in the con- 



tinuum limit. In particular, no other lower- dimensional operators except the one in 
Eq. ( p. 38 ) are to be subtracted non-perturbatively, which is in striking contrast to 



the LQCD formulation with Wilson fermions. Note that compared to other possible 
implementations with staggered fermions, this present method has an advantage of 
the possibility of subtraction at each timeslice. 

Throughout the simulation we use only degenerate mesons, i.e. m s = rrid = m u . 
Since only the negative parity part of O su b contributes in Eq. ( |3.40| ), one naively 
expects infinity when calculating on. However, the matrix elements (0\Oi\K°) of all 
basis operators vanish when m s = rrid due to invariance of both the Lagrangian and 
all the operators in question under the CPS symmetry, which is defined as the CP 
symmetry combined with interchange of s and d quarks [16| . Thus calculation of 



requires taking the first derivative of (0\Oi\K°) with respect to (m^ — m s ). In order to 
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evaluate the first derivative numerically, we insert another fermion matrix inversion in 



turn into all propagators involving the strange quark, as explained below in Sec.|Q 
3.2.3 Contractions To Be Computed 



According to Eqs. (|3.39|) and (|3.40| ), we need to compute three diagrams involving 



four-fermion operators (shown in Fig. |3.2|) and a couple of bilinear contractions. The 
exact expressions for these contractions are given in the following section. Here I 
would like to note that the "eight" contraction type (Fig. |3.2| a) is relatively cheap to 
compute. It is the only contraction needed for the AI = 3/2 amplitude. The "eye" 
and "annihilation" diagrams (Fig. |3.2| b and |3.2| c) are much more expensive since they 



involve calculation of propagators from every point in space-time. Their calculation 
is essential for estimation of the AI = 1/2 amplitude and e'/e . 

3.3 Numerical Techniques 



In Sec. |3.1.5| I have briefly overviewed the conceptual basics of staggered fermions. 
In this section I describe the multitude of methods and tricks used for the purpose of 
calculating the matrix elements of interest in the language of staggered fermions. 

3.3.1 Quark Operators 

There are a number of issues about quark operators that have to be addressed. The 
subject is complicated since each operator is characterized by a given physical flavour, 
staggered spin-flavour and color flow. (Physical flavors should be distinguished from 
staggered flavors; the latter are 4 artificial independent species per each physical 
flavour) . 
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Figure 3.2: Five diagrams types needed to be computed: (a) "Eight"; (b) "Eye"; (c) 
"Annihilation"; (d) "Subtraction"; (e) two-point function. 
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Firstly, there is a choice in constructing four-fermion operators regarding the ar- 
rangement of physical flavour (such as s, d and u) flow. It is related to the fact that 
there are always two ways to arrange bilinears of any given operator, which are re- 
lated by the Fierz transformation. We choose to use such form of the operators that 
all contractions can be represented as a product of two flavour traces; in other words, 
that insertion of the vacuum between the two bilinears would give non-vanishing re- 
sults. To be more precise, for "eight" contractions the operators are rendered in the 
form (s r u)(u T d), while for the "eye" and "annihilation" contractions the appro- 
priate form is (s T d)(q T q). This is done in the continuum, before considering the 
staggered fermion transcription. One should be careful in not forgetting that some 
operators can have multiple non- vanishing contractions of a given type, in which case 
all of them should be cast in the two flavour trace form. For example, the operator 
{s a 7aj(1 — 75)^0) (s/37/x(l -75)5/3) can have non-vanishing two flavour trace "eye" type 
contractions between kaon and pion states in two forms: the present one as well as 
(s«7At(l — 75)^/3) (s/37/i(l — 75)s a ). The contributions of these two contractions should 
be added when calculating the matrix element. 

Secondly, the flavour structure of the operators should be assigned. In general, 
flavour assignment is arbitrary, but in cases like ours when one deals with matrix 
elements of Goldstone bosons, there is no choice left |37| . As mentioned in Sec. p. 1.5 



the external Goldstone bosons states have spin-flavor structure 75 <8> £5. The flavor 
structure of the operators is defined by requiring non-vanishing of the staggered flavor 
traces, and so it depends on the contraction type: the flavor structure is £5 ® £5 in 
"eights", 1 <E> 1 in "eyes". In "annihilation" contractions the flavor structure is 1 
for the bilinear in the internal quark loop trace and £5 for the one involved in the 
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external trace. The bilinears are transcribed in a similar manner: the Goldstone 
boson two-point function has the flavor £5, while the sd operator in the "subtraction" 
contraction has flavor 1. 

Thirdly, either one or two color traces may be appropriate for a particular con- 
traction with a given operator (see the next subsection for details). In one trace 
contractions (type "F") the color flow is exchanged between the bilinears, while in 
two trace contractions (type "U") the color flow is contained within each bilinear so 
that the contraction is the product of two color traces. In either contraction type, 
when the distance between staggered fermion fields being color-connected is non-zero, 
a gauge connector is inserted in the gauge-invariant fashion. The connector is com- 
puted as the average of products of gauge links along all shortest paths connecting 
the two sites. We also implement tadpole improvement by dividing each link in every 
gauge connector by u = (| Tr[[/p]) 1//4 , where Up is the average plaquette value. 
This is done in order to improve the properties of lattice perturbation theory by non- 
perturbatively removing "tadpole diagrams" JITJ. At the same time, in addition to 



rescaling the links, each bilinear should be multiplied by uq to account for quark field 
renormalization. 

3.3.2 Sources And Contractions 

Here we show how sources and sinks are used to construct the correlators, and give 
explicit expressions for various types of contractions which were described in words 
in the previous subsection. 

We use local U(l) pseudofermion wall sources. That is, we set up a field of U(l) 
phases £ Q (x;t ) (|£a| = 1) fc> r each color and each site at a given timeslice to, which 
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are chosen at random and independently, so they satisfy 

(C(x;toH/3(y;to)>=<W^x, y . (3.41) 

(Boldface characters designate spatial parts of the 4-vector with the same name.) 
We proceed to explain how this setup works in the case of the two-point function 
calculation, with trivial generalization to "eight" and "annihilation" contractions. 

Consider the propagator from a wall source at to = in a given background gauge 
configuration, computed by solving the linear system of equations 

J2(P + m)f y X p{y) = £ a (x; 0)4 4 ,o- (3.42) 

y 

In other words, we need to find (p + m) _1 b for an appropriate source b. Note that 
it is equivalent to computing 

xp(y) = E o)GWy; x > °) > ( 3 - 43 ) 

X 

where G(y; x) is the propagator from 4-point x to 4-point y. For staggered fermions 
description we label the fields by hypercube index h and the hypercube corner indices 



A^ e {0, l} 4 , as mentioned in Sec. p. 1.5 



The two-point function is constructed as follows: 

TP = 17 E Xl(h, A)U aP {h, A,A + A)xp(h, A + A) 0(A) (-1) A , (3.44) 

10 h,A 

where phases <f){A) = jTr[T A T s ^a+a anc l distances A M = 2 + are defined 
for a given bilinear operator with spin-flavour quantum numbers S®F; U(h, A, A+A) 
is the appropriate gauge connector (see below), modulo 2 summation is implied for 
hypercube indices A, and h runs over all hypercubes in a given timeslice t where the 
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operator is inserted. The factor (— 1) A takes into account that for staggered fermions 
G(x;y) = G\y; x){— l) x (— l) y - Equation ( |3.44| ) corresponds to 



TP oc £ G a p(z,y) T G^(y,x) (-1)* £(z)£ y (x) , (3.45) 



where T is used for simplicity to show the appropriate operator structure. The sum- 
mation over y is done over the entire space-time. The summation over x and z over 
the spatial volume at the timeslice t averages over the noise, so the last equation is 
equivalent to 

TP oc £ Tr [G(x, y) T G(y, x) {-If} (3.46) 
(compare to Eq. (ft.l5Q ). Therefore, using the pseudofermion wall source is equivalent 



to summation of contractions obtained with independent local delta-function sources. 
Note that the factor (— l) x and zero distance in the staggered fermions language are 
equivalent to spin-flavor structure 75 £§£5. This means the source creates pseudoscalar 
mesons at rest, which includes Goldstone bosons. Strictly speaking, this source also 
creates mesons with spin-flavor structure 7574 ® £ 5 £ 4 , since it is defined only on one 
timeslice instead of two. However, as explained in the first footnote in Section 2.3 of 



Ref. Ij37|| , these states do not contribute. We have used one copy of pseudofermion 
sources per configuration. 

Numerically, the inversion of the (p + m) matrix is done by the conjugate gradient 



algorithm. For details, the interested reader may consult Ref. |42| . 

Analogously to the kaon source, the pion sink at time T is constructed by using 
another set of U(l) random noise ((£*(x;T) ^(y;T)) = 5 a ^ £ Xiy , |£| = 1). The 
propagator $ is computed by solving 

J2(P + m)f v Mv) = T)8 XA}T . (3.47) 
y 
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Suppose Ai, A 2 G {0, l} 4 and (/>i(A), 2 (^4) are distances and phases of the two 
staggered fermion bilinears making up a given four-fermion operator. The expression 
for the "eight" contraction (Fig. |3.2| a) with two color traces ( "U" type) is given by 



E u = -L E Xa(h, A)U a p(h, A,A + A 1 ) X p{h i A + A,) MA) (-1 

iD h,A,B 



x <S>* p (h, B)U pa (h, B,B + A 2 )® a (h, B + A 2 ) fc(B) (-1) B .(3.48) 

In this expression A, B G {0, l} 4 run over 16 corners of a given hypercube (modulo 
2 summation is implied for these indices). The hypercube index h, as before, runs 
over the entire spatial volume of the timeslice t of the operator insertion. The gauge 
connector U(h, A, B) is the identity matrix when A = B, otherwise it is the average 
of products of gauge links in the given configuration along all shortest paths from A 
to B in a given hypercube h. The expression ( 3.48|) , as well as all other contractions, 



is computed for each background gauge configuration (generated beforehand by one 
of the algorithms of Monte Carlo importance sampling) and is subject to ensemble 
averaging. Whenever several contractions are combined in a single quantity, such as 
a B ratio, we use jackknife to estimate the statistical error. 

The expression for one color trace ( "F" type) contraction is similar: 



E f = T& £ Xa(h,A)U af ,(h,A,B + A 2 )xAh,A + A 1 )MA)(-l) 

iD h,A,B 



x$*Jh, B)U p(7 (h, B,A + A^^h, B + A 2 ) fa(B) (-1) B , (3.49) 



For "eye" and "subtraction" diagrams (Fig. |3.2| b and |3.2| d) the source setup is a 



little more involved, since the kaon and pion are directly connected by a propagator. 
In order to construct a wall source we need to compute the product 

1>(y) = X) G(y, t- x, T) ■ G(x, T; x , t )(-l) x '- 

X 
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In order to avoid computing propagators from every point x at the timeslice T, we 
first compute propagator G(x, T; x , to), cut out the timeslice T and use it as the 
source for calculating the propagator to (y,t). This amounts to inverting the system 

Ytf> + m)f y Mv) = Xa(x) <W)(-ir , (3.50) 



where Xa( x ) is the propagator from the wall source at t defined in Eq. ( |3.42| ). We 
use the following expression for evaluating the "subtraction" diagram: 

S = -kj2x* a ( h > A)U a ?(h, A,A + A)Mh, A + A) 0(A) (-1) A , (3.51) 

iD h,A 

Again, averaging over the noise leaves only expressions local in both source and sink, 
so in the continuum language we obtain: 

S oc ]T tr G(x, 0; z, t) T G(z, t; y, T) G(y, T; x, 0) (-l)*(-l) v . (3.52) 

In this work we are mostly interested in subtracting the operator s(l g) l)d, so in 
Eq. (p3ID A = (0, 0, 0, 0) and <p(A) = 1. 



In order to efficiently compute fermion loops for "eye" and "annihilation" diagrams 
(Fig. |3.2jb and |3~2]c), we use U(l) noise copies C , i — 1> • • ',N, at every point in 



space-time. We compute 77 W by inverting (p + m)r]( % > = C . If is easy to convince 
oneself that the propagator from y to x equals 

G{x ] y) = (r ]x Q. (3.53) 

The efficiency of this method is crucial for obtaining good statistical precision: in 
order to obtain accurate results iV should be sufficiently large. In practice we average 
over iV = 10 noise samples per site per color. Note that, since our lattice is copied 2 
or 4 times in the time direction (see next Subsection), independent noise samples are 
used for each of the lattice copies to increase the efficiency. 
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The expression for "U" and "F" type "eye" diagrams are as follows: 
l u = W E X* a (h, A)U aP (h, A,A + AJMh, A + A x ) ^(A) (-1) A 

h,A,B 
1 ^ 

x 77 E C^(h,B)U pa (h,B,B + A 2 )^\h,B + A 2 )MB) (-l) fl (3.54) 
iV i=i 

^ = W E ^^(/i, A, B + A 2 )^(/i, A + Ax) <MA) (-1) A 

1 * 

x 77 E C^(^5)^(^^^ + Ai)^ ) (^5 + A 2 )0 2 (5) (-1) B (3.55) 

JV i=l 

The computation of "annihilation" diagrams (Fig. |3.2| c) is similar to the two-point 
function, except the fermion loop is added and the derivative with respect to the 
quark mass difference m^ — m s is inserted in turn in every strange quark propagator. 
Derivatives of the propagators with respect to quark mass can be obtained by inverting 
equations 

(P + m) X ' = ~X, (3.56) 
(p + m)7]' (i) = -rj® . (3.57) 

Indeed, 

dx_ = _d 1_ = 1 = l_ 

dm dm p + m (p + m) 2 p + m 

We have, therefore, four kinds of "annihilation" contractions. 

Ait, =i^E x:(h,A)U aP (h,A,A + A 1 ) X p(h,A + A 1 ) ( j )l (A)(-l) A 

h,A,B 
1 N 

x 77 E (®*(h,B)U pa (h,B,B + A 2 )r)®(h,B + A 2 )MB) (-1)^3.58) 

- /V i=l 

Aif = lii E Xa(h, A)U aa (h, A,B + A 2 ) X p{h, A + A x ) <pi(A) (-1) A 

h,A,B 
1 N 

x 77 E $ ) *(h,B)U pP (h,B,A + A 1 )r,V(h,B + A 2 )<f> 2 (B) (-1) B (3.59) 

A 2C , =i^E X* a {KA)U aP {h,A,A + A 1 )x^KA + A 1 )<j )l {A){-l) A 

h,A,B 
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1 N 

x mY, C^{h,B)U pa (h,B,B + A 2 )r,'^{h,B + A 2 )<p 2 {B) (-1)^3.60) 



N ■ , 



A 2F = wH X* a (h, A)U aa (h, A,B + A 2 ) X p(h, A + A x ) ^(A) (-1)' 



h,A,B 
N 



1 

X Ci°*C», 5)^/1,5, A + A 1 )iyW(/i,B + A a )0 !1 (S) (-1)^3.61) 



In this subsection I have given explicit expressions for calculating several types of 
contractions. These contractions should be combined in an appropriate way for each 
given operator. This is spelled out in the Appendix. 

3.3.3 Simulation Parameters 



The parameters of simulation are listed in the Tables |3.1| and [3.2| . We use periodic 
boundary conditions in both space and time. Our main "reference" ensemble is a 
set of quenched configurations at (3 = 6/g 2 = 6.0 (Qi). In addition, we use an 
ensemble with the same (3 but a larger volume (Q 2 ), an ensemble with (3 = 6.2 
(Q 3 ) for checking the lattice spacing dependence, and an ensemble with 2 dynamical 
flavors (m = 0.01) generated by the Columbia group, used for checking the impact of 
quenching. In addition, a number of quenched ensembles (Table |3.2| ) with various (3 
are used for taking the continuum limit of Bk and fx- The quenched ensembles were 
obtained using 4-to-l ratio of 3-hit SU(2) overrelaxed and heatbath algorithms. The 
configurations were separated by 1000 sweeps. The dynamical configurations were 
obtained by the R-algorithm [f£|. 

We use renormalized coupling constant which is found as follows. We calculate 
the Lepage-Mackenzie coupling ay from 

1 47T 

- In(-TrUp) = — a v (3.41/a) [1 - (1.185 - 0.070^)^] , (3.62) 
3 3 
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Ensemble 


N f 


P 




Size 






a 1 , 


L, 


Number of 


Quark masses 


name 












GeV 


fm 


conf. 


considered 


Qi 





6.0 


16 3 


x (32 


X 


4) 


2.07 


1.5 


216 


0.01 — 0.05 


Q2 





6.0 


32 3 


x (64 


X 


2) 


2.07 


3.1 


26 


0.01 — 0.05 


Q 3 





6.2 


24 3 


x (48 


X 


4) 


2.77 


1.7 


93 


0.005 — 0.030 


D 


2 


5.7 


16 3 


x (32 


X 


4) 


2.0 


1.6 


83 


0.01 — 0.05 



Table 3.1: Parameters of ensembles used for calculation of (7T7r|Oi_io|i^). 



N f 


p 




Size 






a" 1 , GeV 


L, fm 


Number of 


Quark 


masses 


















configurations 


considered 


2 


5.7 


16 3 


x (32 


X 


4) 


2.0 


1.6 


83 


0.01 - 


- 0.05 





5.7 


16 3 


x (32 


X 


4) 


1.18 


2.7 


259 


0.01 - 


- 0.08 





5.8 


16 3 


x (32 


X 


4) 


1.48 


2.2 


200 


0.01 - 


- 0.04 





5.9 


16 3 


x (32 


X 


4) 


1.77 


1.8 


200 


0.01 - 


- 0.04 





6.0 


16 3 


x (32 


X 


4) 


2.07 


1.5 


221 


0.01 - 


- 0.04 





6.05 


16 3 


x (32 


X 


4) 


2.26 


1.4 


306 


0.01 - 


- 0.05 





6.1 


24 3 


x (48 


X 


4) 


2.44 


2.0 


100 


0.01 - 


- 0.04 





6.2 


24 3 


x (48 


X 


4) 


2.77 


1.7 


121 


0.005 - 


- 0.035 





6.4 


32 3 


x (64 


X 


4) 


3.64 


1.8 


50 


0.005 - 


- 0.030 



Table 3.2: Parameters of ensembles used for calculation of Bx and fx- 



where Nf is the number of sea quark flavors. The coupling in MS scheme can be 
obtained as 

3.41 e 5//6 2 

aj^(3.41/a) = av( — ) [1 H — «v] • (3.63) 

a 7r 

Then using the two- loop running we can obtain at any scale q* . 



As mentioned in Sec. |3.1| , the lattice scale a can be determined by comparing a 
suitable quantity computed on the lattice with the experimental value. In this work 
we use the mass of the p meson for such purposes. 



46 



The quenched lattice scale was set as in Ref. [53], i.e. we demand perturbative 
scaling of the form 



V y MS 7 y MS 7 

normalizing so that the world data for the p mass ^3] is well fit by 

m p (a) = (770 MeV ) ■ (1 + AV(/3)) . (3.65) 

This procedure takes into account the fact that the p mass itself receives corrections 
at finite a. Thus cutoff effects are reduced at any given a. 

The lattice spacing for the dynamical ensemble was also set by the p mass [[46] . 

3.3.4 B Ratios 

It has become conventional to express the results for matrix elements in terms 
of so-called B ratios, which are the ratios of desired four-fermion matrix elements to 
their values obtained by vacuum saturation approximation (VSA). VSA is a crude 
method for obtaining matrix elements of four-fermion operators based on inserting a 
complete set of states between the bilinears and keeping only the vacuum terms. For 
example, the kaon matrix element of the AS = 2 operator can be obtained as 

(1^1(^(1-75^) (s lfl (l- l5 )d) \K°) = ^(K0\(sY(l-l 5 )d)\0) (0|(*y„(l-7fc)d) \K°) 

VSA provides a convenient scale, so that the answers can be quoted as the VSA result 
times a dimensionless quantity called B ratio. 

Sometimes it is convenient to compute B ratios directly on the lattice. For ex- 
ample, the B ratios of operators O2 and O4 are formed by dividing the full ma- 
trix element by the product of axial-current two-point functions (Fig. |3.3j ). We 
expect the denominator to form a plateau in the middle of the lattice, equal to 
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Ze~ mitT (7r|A M |0) • (0| where A^ are the axial vector currents with appropriate 
flavor quantum numbers for kaon and pion. The factor Ze^ 171 ^ cancels, leaving the 
desirable ratio (^10^) / ((7r|A M |0) • (^A^K)). Apart from common normalization 
factors, a number of systematic errors also tend to cancel in this ratio, including 
the uncertainty in the lattice spacing, quenching and in some cases perturbative cor- 
rection uncertainty. Therefore, it is sometimes reasonable to give lattice answers in 
terms of the B ratios. 




Figure 3.3: B ratio is formed by dividing the four-fermion matrix element by the 
product of two-point functions, typically involving A^ or P bilinears. All the operators 
involved are inserted at the same timeslice t, and the external meson sources are 
also located at the same timeslices for numerator and denominator. This enables 
cancellation of various common factors. 

However, eventually the physical matrix element needs to be reconstructed by 
using the known experimental parameters (namely f K ) to compute VSA. In some 
cases, for example in calculation of Bk (see below) , this is easy. In some other cases, 
such as for operators 5 — 8 , the VSA itself is given in terms of (0\P\K), which 
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is known very imprecisely due to the failure of perturbative matching (see Sec. 5.1). 
Then it is more reasonable to give answers in terms of matrix elements in physical 
units. We have adopted the strategy of expressing all matrix elements in units of 
(7r|-(4 M |0) (0|A M |i^) = (f l K tt ) 2m2 M a t an intermediate stage, and using pre-computed 
f l £ tt at the given meson mass to convert to physical units. This method is sensitive 
to the choice of the lattice spacing. 
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Figure 3.4: An example of the signal we get for one of the B ratios (in this case, for 
the "eye" part of the O2 operator on Qi ensemble, using quark mass 0.01). The wall 
sources are at t = 1 and t = 49. We see that the excited states quickly disappear and 
a stable, well-distinguished plateau is observed. We perform jackknife averaging in 
the range of t from 12 to 37 (shown with the horizontal lines). 



It is very important to check that the time distance between the kaon and pion 
sources T is large enough so that the excited states do not contribute. That is, the 
plateau in the middle of the lattice should be sufficiently flat, and the B ratios should 
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not depend on T. We have found that in order to satisfy this requirement the lattice 
has to be artificially extended in time direction by copying the gauge links (4 times in 
the case of the smaller volume lattices, 2 times otherwise). Thus we get rid of excited 
states contamination and wrap-around effects. 
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Figure 3.5: An example of the signal for (0\O2\K°) / [(m^ — m s ) (0|s7s<i|if )] on Qi 
ensemble with quark mass 0.01. The kaon source is at t — 1. We average over the 
range of t from 5 to 12 (shown with horizontal lines). 



We are using T = 72 for Q3 (j3 = 6.2) ensemble, and T = 48 for the rest of the 
ensembles. An example of a plateau that we obtain with this choice of T is shown in 



Figs. [T4] a nd p75| . To read off the result, we average over the whole extension of the 
plateau, and use jackknife to estimate the statistical error in this average. 
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3.4 Computational Summary 

The bulk of the calculations in this thesis consisted of computing fermion propaga- 
tors on a number of gauge configurations and combining the results in ways described 
above in Sec. |3.3|. The first task by far consumes most of computer resources, al- 



though the second one is no less consuming of human labour. Calculations were done 
on high-performance parallel computers CRAY-T3E at Ohio Supercomputing Center 
(OSC) and NERSC. The architecture of CRAY-T3E is very fitting for lattice compu- 
tations since it allows to assign a subregion of the lattice to each processor and run 
the program in parallel, with some exchange of information between the processors 
when necessary. A parallel FORTRAN code was developed by a collaboration of Dr. 
G. Kilcup, Dr. L. Venkataraman and the author. The basic task of inverting the 
fermion matrix (Eq. ( |3.42| )) is done by the conjugate gradient algorithm, adapted to 



the staggered fermion case (discussed in detail in the Ph.D. thesis of L. Venkatara- 
man ^2|). The computer at OSC is capable of providing a peak performance of 600 
MFLOPS (millions of floating point operations per second) per processor. Up to 128 
processors have been used. As is well known, the performance sustained in realistic 
computations is significantly less that the theoretical peak. Our program has been 
able to sustain performance of 100 MFLOPS per processor. The total running time 
taken to collect the data in this thesis is approximately one year. 

In addition to the large-scale calculations on the supercomputer, the raw data 
obtained thereby were analyzed and combined into meaningful results on a UNIX 
workstation. This included statistical analysis, linear and non-linear fitting, solving 
ODE etc. 
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CHAPTER 4 



AI = 1/2 RULE ON THE LATTICE 



Using the data obtained for matrix elements of basis operators, in this section I 
report numerical results for ReA and ReA 2 amplitudes as well as their ratio. These 
amplitudes are discussed separately since the statistics for ReA 2 is much better and 
the continuum limit extrapolation is more reliable. I also present our results for B K 
and f K . 



4.1 AI = 3/2 Amplitude 



The expression for ReA 2 can be written as 



ReA 2 = -^V ud V* s z + (n){0 2 (n)) 



2- 



(4.1) 



where z + (fj) is a Wilson coefficient (we use ji — 2 GeV) and 



(0 2 ) 2 = ((nn) I=2 \0 2 '>\K). 



(4.2) 



Here 




0! 2) = i[(^(l-7 5 )«)(H7 M (l-7 5 M) 

+(*y„(l - 75)^(^(1 " 76)«) - (^(1 - ls)d)(dr(l ~ 7fe)d)]- (4-3) 



52 



In lowest-order chiral perturbation theory the matrix element (02)2 can be expressed 



as 



(0 2 ) 2 = V2 — 2 , (4.4) 

J m M 



so that 



ReA 2 = G F V ud V* s - * R 2 , (4.5) 



where R 2 is calculated on the lattice and is defined as 



8,-™. («) 



m M 

It can be shown that this quantity involves only "eight" diagrams and there are 

(2) 

no subtractions to be made. This follows from the fact that 2 operator belongs 
to the (27,1) representation of the SU{3)l x SU(3)r symmetry group, while any 
contractions of the type "eye" or "subtraction" can give rise only to a term belonging 
to the (8,1) representation. Moreover, in the limit of exact SU (3) navm symmetry R 2 
is directly related E71 to parameter Bk (which is the B ratio of the neutral kaon 



mixing operator, see next subsection), so that 

R 2 =~z + (fi)B K (ii)f 2 K . (4.7) 
4.1.1 Calculating Bk 

The parameter Bk is defined as follows: 

3 m(s lL d) (s~f L d)\K°) 3 (KB\(Tn L d) (Tr L d)\K°) 



(K°\(s lL d)\0) ■ (0\(s lL d)\K°) 8 m K l 



K 



Its calculation is relatively straightforward, based on the techniques presented in 
Chapter |3] (detailed descriptions of the method are also readily available in the lit- 



erature, for example in [25]). It is also relatively cheap to compute since, as already 
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Figure 4.1: Parameter Bk in NDR MS scheme at 2 GeV on the ensemble D vs. 
the meson mass squared. The fit is of the form Bk = a + bm 2 K + c m 2 K \ogm 2 K . The 
vertical line here and in the other plots below marks the physical kaon mass. 



mentioned, it involves only "eight" diagrams". We have computed Bk on all en- 
sembles listed in Table ft.2| . For each ensemble, we study Bk at a variety of quark 
masses, and find the mass dependence well-fitted by the form predicted by the chiral 
perturbation theory 



B K = a + bm 2 K + cm 2 K logm 2 K 
(for example, see Fig. [O for D ensemble). Here we are mostly interested in the value 



of Bk at kaon mass, so there is almost no chiral extrapolation involved. 

Sharpe [^] showed that calculation of Bk with staggered fermions introduces 



scaling violations at worst at the order of a . Our data confirm this prediction (see 
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Figure 4.2: Parameter Bk in NDR MS" scheme at 2 GeV in quenched QCD vs. a 2 . 
Extrapolation to continuum limit is linear, according to Sharpe's prediction. 



Fig. f4.2| ). Extrapolating to the continuum limit, we obtain in quenched QCD: 

B K (NDR, /i = 2 GeV) = 0.554 ± 0.009. (4.9) 

We also found the effect of quenching to be very small, not more than a few percent. 

4.1.2 Calculating fx 

The pseudoscalar decay constant is obtained from the amplitude Ca multiplying 
the exponential in the expression for meson correlator with the axial current operator 
(see Eq. ( p. 1 1| ) ) . The absolute normalization is given by the pion field renormalization 
factor Cp obtained in the analogous manner from the pion wall correlator with the 
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pseudoscalar density operator. Taking various factors into account, one obtains pj|j: 



y/ sinh rriK 



K 



m K /2 



2 sinh J N f V 



(4.10) 




Figure 4.3: Pseudoscalar decay constant on the dynamical (3 = 5.7 (squares) and 
quenched (3 = 6.0 (diamonds) ensembles vs. meson mass squared. The experimental 
number is shown with the horizontal line. 



Fig. fO| shows our results for the pseudoscalar decay constant plotted against the 
mass of the meson squared. Extrapolation to the chiral limit produces the value for 
Thus computing /„- across all quenched ensembles, the continuum limit value of 
134.5 ± 3.6 MeV can be obtained (see Fig. ^4.4| ). Note that this value is allowed to be 
different from the experiment since we consider QCD in the quenched approximation. 
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Figure 4.4: Continuum limit of f n in quenched QCD. 



4.1.3 ReAo Results 



Using our data for B^ and fx, we found that the ratio R 2 shows a large depen- 



dence on the meson mass used in the simulation (Fig. [4.5|) . This is not surprising since 
both Bk and Jk depend on this mass quite significantly. Which meson mass should 
be used to read off the R 2 value for estimation of KeA 2 becomes an open question. If 
known, the higher order chiral terms would remove this ambiguity. Forced to make a 
choice, we extrapolate to M 2 = {m 2 K + m 2 r )/2. 

Taking the continuum limit we obtain in quenched QCD: 

ReA 2 = (1.77 ± 0.07) • 10" 8 GeV, 

where the error is only statistical, to be compared with the experimental result 
ReA 2 = 1.23 • 10- 8 GeV. 
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Figure 4.5: Matrix element R2 computed on the dynamical D (squares) and quenched 
Qi (diamonds) ensembles. ReA 2 is proportional to this quantity at the lowest order 
in chiral perturbation theory. The horizontal line shows the value corresponding to 
the experiment. 
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Higher order chiral terms (including the meson mass dependence) are the largest 
systematic error in this determination. According to Golterman and Leung EDI , one- 
loop corrections in (quenched) chiral perturbation theory are expected to be as large 
as 30 — 60%. Other uncertainties (from lattice scale determination, from perturbative 
operator matching and from using finite lattice volume) are much smaller. 

4.2 AI = 1/2 Amplitude 

Using Eqs. ( |3.39| , |3.40"D , ReA can be expressed as 

ReA = ^V ud V: s ^^R , (4.11) 

where 



-n-0 = 2^ Zi 2 



Here Z{ are Wilson coefficients, and the subscript V indicates that these matrix 
elements already include subtraction of (tt + \O su i ) \K + } . Of^ are isospin parts of 
operators 0{ (given in the Appendix |A| for completeness). For example, 

0? ] = ^(^(l-7fe)^(^(l-7b)«)-^(l-7fe)«)(^(l-7fe)^ 

+ l(s lft (l- l5 )d)(d^(l- l5 )d) (4.12) 

{ ? = ^(s7, t (l-75)«)(«7 M (l-75)rf)-^(^(l-75)rf)(^(l-75)«) 



+ ~(s Tm (1 - l5 )d)(dY(l - 7s)d) (4-13) 



All contraction types are needed in order to calculate Ro (as opposed to the case of 
calculation of ReA 2 ), including the expensive "eyes" and "annihilations". The results 
for quenched (3 = 6.0 and (3 = 6.2 ensembles are shown in Fig. [4.6| . Dependence of Ro 
on the meson mass is small, so there is no big ambiguity about the mass prescription 
as in the R 2 case. 
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Figure 4.6: Matrix element Rq for quenched ensembles plotted against the meson 
mass squared. ReA is proportional to this quantity in the lowest order in chiral 
perturbation theory. The upper group of points is for ensembles Qi and Q 2 , while 
the lower group is for Q 3 . Only statistical errors are shown. The horizontal line shows 
the value corresponding to the experiment. 
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Figure 4.7: ReA for quenched ensembles plotted against lattice spacing squared. A 
naive extrapolation to the continuum limit is made. The horizontal line corresponds 
to the experimental result of 27.4 • 10~ 8 GeV. Only statistical errors are shown. 



It is apparent from Fig. [4.7| that considerable cutoff dependence is present. The 
continuum extrapolation shown in the figure is, naturally, not as reliable as one would 
desire. In the future, data can be collected on a few more ensembles to allow a more 
reliable extrapolation. Note that I have plotted KeA vs. a 2 . Strictly speaking, it is 
not known whether the leading order cutoff effects are 0(a) or 0(a 2 ). For the case 
of Bk it was shown |I9| that these effects are 0(a 2 ), and we have assumed it is still 
the case for ReAo- In any case, with present data the precision of the extrapolation 
is far from the level where the form of the fit matters significantly. 
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The effect of the final state interactions (contained in the higher order terms of 
the chiral perturbation theory) is likely to be large. This is the biggest and most 
poorly estimated uncertainty. 
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Figure 4.8: Comparison of quenched (Qi) and dynamical results for ReA at compa- 
rable lattice spacings. 



In addition, there is an operator matching uncertainty coming from mixing of O2 
with Oq operator through penguin diagrams in lattice perturbation theory. This is 



explained in the Section [54]. We estimate that this uncertainty does not exceed 20% 
for all ensembles. 

As for other uncertainties, we have checked the lattice volume dependence by 
comparing ensembles Q\ and Q2 (1.6 and 3.2 fm at (3 — 6.0). The dependence was 
found to be small, so we consider (1.6 fm) 3 as a volume large enough to contain 
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the system. We have also checked the effect of quenching and found it to be small 
compared to noise (see Fig. |4~8| ). 
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Figure 4.9: Contribution of various operators to KeA^. 

The breakdown of contributions of various basis operators to KeA is shown in 
Fig. [4.9| . By far, 2 plays the most important role, whereas penguins have only a 
small influence. 

4.3 Amplitude Ratio 



Shown in Fig. |4.10| is the ratio Re^/ReAo as directly computed on the lattice 
for quenched and dynamical data sets. The data exhibit strong dependence on the 
meson mass, primarily due to the chiral behaviour of KeA 2 (compare with Fig. |4.5| ). 
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Figure 4.10: ReA 2 /ReA versus the meson mass squared for quenched and dynamical 
ensembles. Ensembles Q± and D have comparable lattice spacings. The dynamical 
ensemble data were used for the fit. The horizontal line shows the experimental value 
of 1/22. The error bars show only the statistical errors obtained by jackknife. 
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4.4 Summary 

Within our errors the results seem to confirm, indeed, the common belief that most 
of the AI = 1/2 enhancement comes from strong interactions, in particular from the 
"eye" and "annihilation" diagrams. The exact amount of enhancement is broadly 
consistent with experiment while being subject to considerable uncertainty due to 
higher-order chiral terms. Other systematic errors are the same as those described in 



Section 4.2 
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CHAPTER 5 



OPERATOR MATCHING AND e'/e 



As mentioned before, we have computed the matrix elements of all relevant op- 
erators with an acceptable statistical accuracy. These are regulated in the lattice 
renormalization scheme. To get physical results, operators need to be matched to 
the same scheme in which the Wilson coefficients were computed in the continuum, 



namely MS* NDR. While perturbative matching works quite well for Rev4 and Re^, 
it seems to break down severely for matching operators relevant for e'/e. 

5.1 Perturbative Matching 

Conventionally, lattice and continuum operators are matched using lattice pertur- 
bation theory: 

OrV) = Of + $>« Hi**) + C i5 )Of + 0(g 4 ) + 0(a n ), (5.1) 

where 7^ is the one-loop anomalous dimension matrix (the same in the continuum 
and on the lattice), and CV,- are finite coefficients calculated in one- loop lattice pertur- 
bation theory ]50, 51]. We use the "horizontal matching" procedure [|52], whereby the 



same coupling constant as in the continuum (g^g) is used. The operators are matched 
at an intermediate scale q* and evolved using the continuum renormalization group 
equations to the reference scale /1, which we take to be 2 GeV. 
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(a) 

(b) 




(d) 



Figure 5.1: Example of four kinds of diagrams with arbitrary number of loops arising 
in renormalization of four-fermion operators: in (a) and (b) no propagator crosses 
the box or the axis; (c) and (d) exemplify the rest of the diagrams. The rectangular 
drawn in dotted line in (b) corresponds to operator structure PPeu- 
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5.1.1 Perturbative Matching And Re^4o,2 

In calculating of KeA Q and KeA 2 , the main contributions come from left-left op- 
erators. One-loop renormalization factors for such (gauge- invariant) operators were 
computed by Ishizuka and Shizawa |50| (for current-current diagrams) and by Patel 
and Sharpe |51] (for penguins). These factors are fairly small, so at the first glance 
the perturbation theory seems to work well, in contrast to the case of left-right oper- 
ators essential for estimating s'/e, as described below. However, even in the case of 
ReA there is a certain ambiguity due to mixing of O2 with 6 through penguin dia- 
grams. The matrix element of Oq is rather large, so it heavily influences (O2) in spite 
of the small mixing coefficient. The operator Oq receives enormous renormalization 
corrections in the first order, as discussed below. Therefore, there is an ambiguity as 
to whether the mixing should be evaluated with renormalized or bare 0§. That is, 



the higher-order diagrams (such as Fig. |5.1| b and |5.1| d) may be quite important here 



In order to estimate the uncertainty of neglecting higher-order diagrams, we eval- 
uate the mixing with 0§ renormalized by the partially non-perturbative procedure de- 
scribed below, and compare with results obtained by evaluating mixing with bare Oq. 
The first method amounts to resummation of those higher-order diagrams belonging 



only to type (b) in Fig. |5.1| , while the second method ignores all higher-than-one- 
order corrections. Results quoted in the previous Section were obtained by the first 
method, which is also close to using tree-level non-diagonal matching. The second 
method would produce values of ReA Q lower by about 20%. Thus we consider 20% a 
reasonable estimate of the matching uncertainty. 

In calculating e'/e the operator matching issue becomes a much more serious 
obstacle as explained below. 
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Figure 5.2: Three contributions to (ti + \Oq\K + ): "eight" (boxes), "eye" (diamonds) 
and "subtraction" (crosses). These data represent bare operators for the dynamical 
ensemble. The fit is done for the sum total of all contributions. 



5.1.2 Problems With Perturbative Matching 

The value of e' / e depends on a number of subtle cancellations between matrix el- 
ements. In particular, Oq and 0§ are considered the most important operators whose 
contributions have opposite signs and almost cancel. Furthermore, the matrix ele- 
ments of individual operators contain three main components ("eights", "eyes", and 



"subtractions"), which again conspire to almost cancel each other out (see Fig. |5.2| ). 
Thus e'/e is extremely sensitive to each of these components, and in particular to 
their matching. 

Consider fermion contractions with operators such as 

(PP)eu = (375®&«)(E75®&d) (5-2) 
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Quark mass 


0.01 


0.02 


0.03 


Bare 


-2.95 ±0.27 


-2.61 ± 0.11 


-2.50 ±0.08 


q* = n/a 


-0.69 ±0.10 


-0.66 ±0.04 


-0.66 ±0.03 


q* = 1/a 


0.01 ±0.07 


-0.07 ±0.03 


-0.11 ±0.02 


Partially non-pert. 


2.03 ±0.07 ±0.42 


0.74 ±0.03 ±0.37 


0.44 ±0.02 ±0.31 


Quark mass 


0.04 


0.05 




Bare 


-2.60 ±0.16 


-2.94 ±0.27 




q* = n/a 


-0.73 ±0.06 


-0.83 ±0.09 




q* = 1/a 


-0.19 ±0.04 


-0.20 ±0.03 




Partially non-pert. 


0.23 ±0.03 ±0.28 


0.25 ±0.05 ±0.36 





Table 5.1: (tt + 7t~\Oq\K°) for Q\ ensemble in units of (GeV) 3 with tree-level matching 
(bare); with one- loop perturbative matching using two values of q*; and with match- 
ing obtained by the partially non-perturbative procedure. The errors are statistical, 
except for the last line, where the second error comes from uncertainty in our deter- 
mination of Zs and Zp. As mentioned in text, there is an unknown uncertainty in 
the partially non-perturbative procedure. 



(SS) W = (sl®ld)(dl®ld) (5.3) 
(PS) A 2u = (375®&d)(dl®ld), (5-4) 

which are main parts of, correspondingly, "eight", "eye" and "subtraction" compo- 
nents of Oq and Os (see the Appendix). The finite renormalization coefficients for 



these operators have been computed in Ref. ||51|| . The diagonal coefficients are very 

large, so the corresponding one-loop corrections are in the neighborhood of —100% 

and strongly depend on which q* is used (refer to Table |5.1| ). Thus perturbation 

theory fails in reliably matching the operators in Eqs. ( |5.2| ) — ( |5.4| ). 

The finite coefficients for other (subdominant) operators, for example (PP)ef, 

(SS)eu and (SS)ef, are not known for formulation with gauge-invariant operators^. 
3 Patel and Sharpe [[5i"| have computed corrections for gauge- noninvariant operators. Operators 



in Eqs. ( |5.2| ) — (5.4) have zero distances, so the corrections are the same for gauge invariant and non- 
invariant operators. Renormalization of other operators (those having non-zero distances) generally 
differs from the gauge-noninvariant operators. 

70 



For illustration purposes, in Table fj]we have used coefficients for gauge non- invariant 



operators computed in Ref. but strictly speaking this is not justified. 

To summarize, perturbative matching does not work and some of the coefficients 
are even poorly known. A solution would be to use a non-perturbative matching 
procedure, such as described in Ref. . Before we proceed to this procedure, let me 
discuss a simpler procedure for non-perturbative calculation of matching coefficients 
for bilinears Z$ and Zp, which is important in its own right since the perturbative 
matching error is the dominant source of errors in determination of light quark masses 
via lattice QCD. 

5.2 Nonperturbative Matching Of The Scalar And Pseudoscalar 
Bilinears And Determination Of Light Quark Masses 

5.2.1 Framework 

The general method of nonperturbative operator matching has been suggested 
by Martinelli et al. ||54j| . It has been successfully applied in several calculations ||53|| . 
Here we follow one variant of the method which allows only the determination of Z$ 
and Zp renormalization factors of the scalar and the pseudoscalar operators, which 
is enough for our purposes (see also Ref JB|)- 

The method is based on averaging the quark propagator between external single 
quark states over gauge configurations in a fixed (usually, Landau) gauge, considering 
off-shell momenta. The momenta should be large enough {\p 2 \ 3> Aq CD ) so that it 
makes sense even to consider single quark states, and on the other hand, the momenta 
should be sufficiently small compared to 1/a in order to avoid large discretization 
effects. Note that these criteria apply to any (perturbative or non-perturbative) 
method of operator matching. 
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The ensemble-averaged propagator is expected to have the form 

{S{P)) = M(p) + izJi%i)MP»a) ' (5 ' 5) 
where Z q is the renormalization of the fermion field. It can be shown that 

Zsip) = (5.6) 

Z P (P) = ^, (5.7) 
m 

where m is the bare quark mass and Zs and Zp are the bilinear renormalization 
factors defined by 

(S, P) latt = Z s ,p (S, P) cont . (5.8) 
It follows that M(p) can be found as 

M W = MM (5 . 9) 

where the trace is taken with respect to both color and Dirac indices. Z q can be 
obtained from 

t = _. o E^ISM-U^i)] ^ (510) 

E M Tr[l] sm (p M a) 

This procedure amounts to resummation of all perturbative corrections to Zs and 
Zp. In the continuum, it corresponds to the so-called regularization-invariant (RI) 
scheme. In order to obtain results in the MS scheme, Zs and Zp need to be rescaled by 



an appropriate factor Z c . This factor was computed at two- loop order in Ref. |p5[ . It 
is of the order of 0.89 for the ensembles of interest, and the higher-order perturbative 
corrections can be safely neglected. Thus Zs t p in MS scheme are given by: 

ZB^) = U us(^ q*)Z c {q*)Zf P {q*), (5.11) 

where Uy^ is the three-loop evolution factor. 
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5.2.2 Procedure 

In applying the general method described above to the case of staggered fermions, 
care needs to be taken in considering the propagator. The complication here is that 
staggered quark fields "live" on the coarse lattice (2a) and there are several different 
but equivalent ways to write down the action. Here I carefully describe the procedure 
and ideas behind our computations. 



Recall that there are two equivalent ways to work with staggered fermions [p26| , [35 
The first one is based on fermion fields x( x : A) defined in the coordinate space (as in 
Section |3.1.5|) . 

Another way is to start with fields ip{k) defined completely in momentum space. 
As is well known, there are 16 poles in the fermion propagator at p^ = (0, n/a). 
These poles are used to identify fermion spin and flavor. Thus it is convenient to 
divide the Brillouin zone into 16 regions, labeled by H G {0, l} 4 . The momentum 
k can be broken down into a physically meaningful part p G [— n/2a, 7r/2a[ 4 and the 
offset corresponding to the given region as follows: k = p + -H. The (free) fermionic 
action can be written as 

S F = 77 E E F )l m *f,h + ~ E sin M (J»®1) F , H MP, F). (5.12) 

V F,H a H 



The matrices (75 ® £f) are related to (75 ® ^f) by a unitary transformation: 

(75 ® t F ) AB = e t\(-) a ' c (75 ® zf) cd h d - b , 

CD 10 

where A ■ C = 5Z„ A^C^. The fields i[)(k) are Fourier transforms of the fields x{ x -> A) 
(i/j(k) = J2 x ,a exp(— ik(x + aA))x(x, A)). Consider fields 4>(p, A) defined by 

<p(p, A) = exp(— ip(x + aA))x(x, A) . (5.13) 

X 
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They are related to if) by 

iP(p,H)=Y,H A - H ^P,A) . (5.14) 

A 

The (free) action in terms of the new fields is given by 

16 _ i 

Sf = v E E 0(P» B ) [ m W + - E sin (P/, a ) (7m ® l) BiA ] 4>(P, A). (5.15) 

^ P A.B a M 



Actions in coordinate and momentum space are completely equivalent |35| , |26| , |56~ 



It is with the latter (momentum space) form of action that we choose to work here. 

In practice, we use 16 5-function sources at the corners of the hypercube contain- 
ing the origin (x = (0, A), A 6 (0, l) 4 ). Propagator G(y, B\ 0, A) is computed for each 
space-time point, where y denotes a hypercube and B is a hypercube corner. We 
perform Fourier transform in y to obtain G(p, B, A), where p G [— 7r/2a, 7r/2a[ 4 . In 
addition, in order to work with fields <p(p, A), G(p,B,A) is multiplied by the phase 
exp [—iap(B — A)] to obtain S(p, B, A). The propagator is averaged over the gauge 
ensemble (denoted () below), and afterwards the inverse is taken. The inverse prop- 
agator then has the form (up to various factors that can be absorbed into Z q ) 

(S^B^A))- 1 = (E ~ sin(ap^) ( 7 „ ® 1) B , A + M{p)Jl®l) BA ) / Z q 



(E - sin(ap M ) r^(A) 5 (A , B+A) + Af(p) 5 A>B ) / Z q . (5.16) 



In 5(A,B+fi) the summation is understood modulo 2. Thus Z 9 can be obtained as 
„-i ■ J2aJ2^(S{p,A,A + fi))- 1 !]^) sm{p^a) 

Z q = ~™ " ^ ■ 27 N > 5 " 17 



and 



M(p)=E(%»A^))~%- (5.1* 
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Here (S(p, A, B))^ 1 is the inverse of the matrix (S(p))ab i n both color and hypercube 
corner index space. In practice, the calculations are made easier by the fact that the 
propagator is proportional to the identity matrix in color space when one works in 
Landau gauge. 

A number of checks have been performed to confirm the expected form of the 
propagator in Eq. ( |5.5|) . In particular, the traces of the averaged inverse propagator 
with all spin-flavor structures other than (7^ ® 1) and (1 <S> 1) were found to be 
consistent with zero within given statistical errors, and the momentum dependence 
of the (7^ ® 1) trace was found to be extremely close to sin (ap^). 

0.20 
0.15 
s 0.10 
0.05 
0.00 

0.00 0.02 0.04 0.06 

ma 

Figure 5.3: M(p) vs. quark mass for one momentum withp 2 = (2.2 GeV) 2 at (5 = 6.0. 
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Figure 5.4: Zp at ma = 0.01 vs. momentum at f3 — 6.4. 
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5.2.3 Results For Z s And Z P 

We have computed Zg and Zp on three quenched and one dynamical ensembles 
(see Table [T2]for details). 



N f 


P 


Size 


or 1 , GeV 


L, fm 


Number of 


Quark masses 










configurations 


considered 


2 


5.7 


16 3 x 32 


2.0 


1.6 


49 


0.02, 0.03 





6.0 


16 3 x 32 


2.07 


1.5 


70 


0.01, 0.02, 0.03, 0.10 





6.2 


24 3 x 48 


2.77 


1.7 


69 


0.01, 0.015, 0.02, 0.03 





6.4 


32 3 x 64 


3.64 


1.8 


57 


0.0100, 0.0125, 0.0150 



Table 5.2: Parameters of ensembles used for calculation of Zs, Zp and light quark 
masses. 



Plotted in Fig. IOI is M(jp) at a given p with q* = \fp 2 = 2.2 GeV, as a function of 



quark mass at /3 — 6.0. The dependence is linear, which is also the case for the rest of 
the momenta. From Eq. (|5.6| ) it is clear that Zs is the slope of this function. The fact 
that the intersect is non-zero signals the presence of non-perturbative pseudoscalar 
vertex contribution, arising from the capability of the pseudoscalar density operator 
to create Goldstone bosons from the vacuum. This contribution is a consequence of 
spontaneous chiral symmetry breakdown, and corresponds to the dynamical genera- 
tion of a (constituent) light quark mass, of magnitude about 300-400 MeV, off-shell 



and gauge-dependent (for further discussion see Refs. [19|, [57| and references therein). 
This effect (measured by the intersect of M(p) as a function of m) behaves as ~ ^ 
for large enough p 2 (p 2 ^> m 2 ). 

Since the intersect is non-zero, Zp is not identical to Zg. This is in contrast to 
Zg = Zp in perturbation theory, where Goldstone boson pole is entirely omitted. 
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The ratio of Zp/Zg (Fig. |5.5|) is almost constant and close to 1 for larger momenta, 
while it grows rapidly at smaller momenta. This shows that the Goldstone boson 
pole contribution is very large. It is proportional to 1/m (again, see Eq. ( |5.7| )). Thus 
Zp (Fig. diverges in the limit of small momenta and quark masses. 

Zs is finite and independent of quark mass (Fig. |5.6|). The errors were obtained 



by combining line fitting and jackknife. In Fig. [5J| I compare the non-perturbative 
results for Zs with those obtained in improved perturbation theory (using gj^g). 

We have checked the discrete rotation invariance by comparing results for several 
momenta related by rotation by 90° around one or several of the axes. Another 
interesting check is to compare results obtained with two momenta which have the 
same p 2 but different component structure, for example: (1,1,1,1) and (0,0,0,2). This 
allows one to check for continuous rotation invariance. As expected, this invariance 
is slightly broken by the lattice cutoff. In the above figures, all kinds of momenta 
component structures are present. Sometimes, two (or more) momenta have the same 
p 2 . By looking at the differences in the plotted function between these two momenta 
it is possible to obtain an indication of the discretization effects. Indeed, these effects 
grow larger for larger momenta, as expected. In addition, they decrease for any given 
momentum as f3 grows (i.e., closer to the continuum limit). These effects are discussed 
more in the next subsection. 

5.2.4 Results For Light Quark Masses 

Values of light quark masses can be obtained from lattice calculations with stag- 
gered fermions as follows. By computing masses of mesons at a variety of bare quark 
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Figure 5.5: Ratio Z P /Z S for (5 = 6.4 vs. momentum. 
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Figure 5.6: Z s at f3 — 6.4 vs. momentum. Solid curve shows the Z s value calculated 
in perturbation theory with g-^ coupling. 
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mass values, the values of bare quark masses are obtained by extrapolation to the cor- 
responding physical meson masses. For example, the strange quark mass corresponds 
to the lattice meson mass being a/2 times the physical kaon mass (this is because the 
kaon on the lattice is composed of two equal-mass quarks, while in the real world the 
kaon is composed of the strange quark and the comparatively massless down-quark; 
and m 2 K oc (m s +m d )). The bare quark mass in physical units is obtained by dividing 
the mass in lattice units by a: 

m = {ma)M ^ ■ (5.19) 
a 

This procedure gives the masses of the quarks in lattice regularization scheme. But 
their values need to be matched to continuum just the same as for the operators, and 
herein lies a well-known problem, since the perturbative operator matching does not 
work very well in this case. 

To overcome this problem, we use the bilinear renormalization constant Zs, com- 
puted non-perturbatively as described in the previous subsection. Using 

z s = Z~ , 

we compute the quark masses in MS NDR scheme as follows: 

m cont (MS NDR, /i) = U{(i,q*)Z c {q*)Zf l {q*)m liitt . (5.20) 



In Fig. |5]7] I show the results for m s at the scale of 2 GeV in MS NDR scheme 
obtained from a variety of matching momenta at f3 = 6.4. Ideally, one expects to see a 
plateau in the middle of the plot, where p 2 is sufficiently above the hadronization scale, 
and sufficiently low so that p 2 a 2 <C 1. However, the data show (1) sizable dependence 
on the momentum component structure (a spread along the vertical axis) and (2) 
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Figure 5.7: Strange quark mass at 2 GeV {(5 = 6.4) vs. the matching momentum 
scale p 2 = q* 2 . 
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Figure 5.8: Strange quark mass at 2 GeV (f3 = 6.4) vs. the matching momentum 
scale p 2 = q* 2 for the axial and diagonal momentum families. 
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a noticeable downward trend. Both are evidently due to the discretization effects, 
and are expected to be proportional to p 2 a 2 , which is not unreasonable according 



to the data. In order to study it more, in Fig. [5]8] I have plotted two families of 
momenta each having unique component structure but different overall magnitude. 
These families are: (1) axial, with momenta (0,0,0,1), (0,0,0,2), (0,0,0,4) etc.; and 
(2) diagonal, with momenta (1,1,1,1), (2,2,2,2) etc. One expects these families to be 
the opposite extreme cases, with the rest of the momenta being in between them in 
certain sense. Linear extrapolation to p 2 = gives very close results for these two 
families. Therefore, such an extrapolation is the most natural way to obtain reliable 
estimates of the quark mass at any given (3, since a large part of the discretization 
effects is removed. 

One can also study the discretization effects by performing continuum extrapola- 



tion. This is done in Fig. |5.9| for ensembles with (3 = 6.4, 6.2 and 6.0. We are using 
the diagonal family for the extrapolation, since arguably it has minimal discretization 
corrections. Two sets of data are plotted in Fig. |5.9|: one corresponding to using the 



matching scale p 2 = 22 GeV 2 , and the other one obtained at each (3 by extrapolating 
to p 2 = 0, as discussed above. The results are consistent, which was expected since 
the discretization errors vanish in continuum limit. This adds to the strength of the 



hypothesis that the slope and spread in data in Fig. |577] is due to discretization effects. 

As our best estimate we quote the value of the strange quark mass in quenched 
QCD0: 

m s (MS NDR, 2 GeV, a = 0) = 102.9 ± 5.1 MeV. (5.21) 

4 In obtaining this result we have used kaon mass to fix the bare strange quark mass, as opposed 
to another frequently-used alternative of using the 4> meson mass 
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0.0 0.1 0.2 

a 2 , GeV 2 



Figure 5.9: Strange quark mass at 2 GeV, a continuum extrapolation in quenched 
QCD. Continuum limit results are shown with bursts at a 2 = 0. 



This was obtained by extrapolating to p 2 = at each (3. The mass m = (m u + md)/2 
can be obtained by dividing the above value by 26, according to chiral perturbation 



theory. Our results agree with those of the JLQCD group [IE], which were obtained 



with a slightly different nonperturbative method. These results are also consistent 



with those obtained with O(a) improved Wilson fermions [20 1 . To close the discussion 



of light quark masses, I report the result for dynamical Nf = 2 ensemble at f3 = 5.7: 
m s (MS NDR, 2 GeV) = 105.1 ± 16.4 MeV. (5.22) 
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At current level of precision, the dynamical result is consistent with the quenched 
ones. 

5.3 Partially Nonperturbative Matching For Basis Opera- 
tors. 

Ideally, one would like to perform a matching procedure for four-fermion opera- 
tors, not only for bilinears. It is quite conceivable that such a procedure will be done 
in the nearest future. As a temporary solution, however, we have adopted a partially 
nonperturbative operator matching procedure, which makes use of bilinear renormal- 
ization coefficients Z P and Z s computed in the previous section. An estimate of the 
renormalization of four-fermion operators can be obtained as follows. 




(b) 



Figure 5.10: Example of one loop diagrams arising in renormalization of four-fermion 
operators: in type (a) no propagator crosses the axis, and type (b) includes the rest 
of the diagrams. 
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Consider renormalization of the pseudoscalar-pseudoscalar operator in Eq. ( |5.2|) . 
At one-loop level, the diagonal renormalization coefficient Cpp (involving diagrams 
shown in Fig. |5.10| ) is almost equal to twice the pseudoscalar bilinear correction Cp. 
This suggests that, at least at one-loop level, the renormalization of (PP)eu comes 
mostly from diagrams in which no gluon propagator crosses the vertical axis of the 



diagram (for example, diagram (a) in Fig. |5 . 1 0|) , and very little from the rest of the 



diagrams (such as diagram (b) in Fig. |5.10|) . In other words, the renormalization of 



(PP)eu would be identical to the renormalization of product of two pseudoscalar 
bilinears, were it not for the diagrams of type (b), which give a subdominant contri- 
bution. Mathematically, 

(PP) c ^ t = (PP) l E ^Z P p + ..., 

with 



Zrr = Z 2 P (1 + ^Cp P + 0(p 4 )) ; (5.23) 



Z P = l + ^Cp + 0(g 4 ) 1 (5.24) 
and dots indicate mixing with other operators (non-diagonal part). The factor Cpp = 



Cpp — 2Cp contains diagrams of type (b) in Fig. |5.10| and is quite small 



In order to proceed, it may be reasonable to assume that the same holds at 



all orders in perturbation theory; namely the diagrams of type (c) in Fig. |5.1| give 
subdominant contribution compared to type (a) of the same Figure. This assump- 
tion should be verified separately by performing non-perturbative renormalization 
procedure for four-fermion operators. If this ansatz is true, we can substitute the 
non-perturbative value of Zp into Eq. ([5.23 ) instead of using the perturbative ex- 



pression from Eq. ( |5.24| ). Thus a partially nonperturbative estimate of (P?)™ n( is 
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obtained. This procedure is quite similar to the tadpole improvement idea: the bulk 
of diagonal renormalization is calculated non-perturbatively, while the rest is reliably 
computed in perturbation theory. Analogously we obtain diagonal renormalization 
of operators (SS) W and (PS) A 2U by using Z S s = Z 2 s (l + jf^rCss + 0(g 4 )) and 
Zp S = Z S Z P (1 + j^Cps + 0(g 4 )). (Note that Z P ^ Z s , even though they are 
equal in perturbation theory.) We match operators at the scale q* = 1/a and use the 
continuum two-loop anomalous dimension to evolve to /i = 2 GeV. 

Unfortunately, the above procedure does not solve completely the problem of 
operator renormalization, since it deals only with diagonal renormalization of the zero- 



distance operators in Eqs. ( |5.2| ) — ( |5.4|) . Even though these operators are dominant in 
contributing to e'/e, other operators (such as (SS)eu and (PP)ef) can be important 
due to mixing with the dominant ones. The mixing coefficients for these operators 
are not known even in perturbation theory. For a reasonable estimate we use the 
coefficients obtained for gauge non-invariant operator mixing pHl . 



Secondly, since renormalization of operators (PP) EU , (SS)iu and (PS)a2U is dra- 
maticf], their influence on other operators through non-diagonal mixing is ambiguous 
at one-loop order, even if the mixing coefficients are known. The ambiguity is due to 



higher order diagrams (for example, those shown in Fig. |5.1| ). In order to partially 
resum them we use operators (PP)eu, (SS)iu and (PS)a2U multiplied by factors Z P) 
Zg and ZpZs, correspondingly, whenever they appear in non-diagonal mixing with 
other operators. This is equivalent to evaluating the diagrams of type (a) and (b) 



in Fig. |5.1| at all orders, but ignoring the rest of the diagrams (such as diagrams (c) 

and (d) in Fig |5.1| ) at all orders higher than first. A completely analogous procedure 

5 For example, at m q = 0.01 and /i = 2 GeV for Qi ensemble we obtain Zpp = 0.055 ± 0.007, 
Z PS = 0.088 ± 0.007 and Z ss = 0.142 ± 0.010. 



was used for mixing of 0$ with O2 through penguins when evaluating KeA . To 
estimate a possible error in this procedure we compare with a simpler one, whereby 
bare operators are used in non-diagonal corrections (i.e. we apply strictly one-loop 
renormalization) . The difference in e'/e between these two approaches is of the same 
order or even less than the error due to uncertainties in determination of Zp and Z$ 
(see Tables |573| and WM . 



Quark 


0.01 


0.02 


0.03 


mass 








Ml (p.r.) 


-61.2 ±2.8 ±10.6 


-27.4 ±0.9 ±8.9 


-16.8 ±0.5 ±8.0 


Ml (1-loop.) 


-52.3 ±2.2 ±10 


-22.0 ±0.8 ±8.3 


-12.2 ±0.5 ±6.9 


M2 (p.r.) 


-38.6 ±2.1 ±6.0 


-18.7 ±0.3 ±7.0 


-11.7±0.2±6.0 


M2 (1-loop) 


-45.4 ±3.5 ±8.6 


-18.8 ±0.4 ±7.0 


-10.3 ±0.3 ±6.0 


M3 (p.r.) 


-97 ± 14 ± 13 


-81 ±4 ±23 


-79 ± 2 ± 27 


M3 (1-loop) 


-142 ±28 ±29 


-88 ± 5 ± 35 


-75 ± 2 ± 39 


Quark 


0.04 


0.05 




mass 








Ml (p.r.) 


-8.0 ±0.9 ±7.2 - 


-4.4 ±0.9 ±7.2 




Ml (1-loop) 


-4.2 ±1.1 ±6.5 - 


-1.2 ± 1.0 ±6.6 




M2 (p.r.) 


-6.1 ±0.5 ±5.3 - 


-3.1 ±0.5 ±4.9 




M2 (1-loop) 


-3.7 ±0.8 ±5.8 - 


-0.9 ±0.8 ±5.2 




M3 (p.r.) 


-81 ±5±38 


-74 ± 4 ± 38 




M3 (1-loop) 


-64 ± 4 ± 52 


-55±5±51 





Table 5.3: e'/e in units of 10~ 4 for Qi ensemble {(3 = 6.0), computed in three ways 
mentioned in text. In all cases, partially- nonperturbative matching have been used 
to obtain the results. 



5.4 Estimates Of e'/e 

Within the procedure outlined in the previous section we have found that (0 6 ) 
has a different sign from the expected one due to a large renormalization factor. This 
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Quark mass 


0.005 


0.010 


0.015 


Ml (p.r.) 


-68.1 ±6.9 ±36.0 


-33.6 ±2.9 ±23.9 


-24.9 ± 1.8 ±22.0 


Ml (1-loop) 


-60.9 ±6.9 ±31.2 


-29.3 ±2.9 ±21.0 


21.4 ± 1.9 ± 19.4 


M2 (p.r.) 


-43.9 ±9.5 ±16.5 


-33.1 ±3.8 ± 18.3 


-22.1 ± 1.2 ± 16.4 


M2 (1-loop) 


-53.6 ±16.9 ±25.0 


-37.9 ±6.0 ±25.3 


-23.0 ± 1.5 ±20.0 


M3 (p.r.) 


-63.3 ±35.1 ± 15.5 


-103.9 ±31.7 ±45.2 


-82.4 ± 12.9 ±51.0 


M3 (1-loop) 


-98.3 ±72.9 ±46.2 


-138.1 ± 53.1 ± 101.5 


-92.4 ± 16.3 ±86.3 


Quark mass 


0.020 


0.030 




Ml (p.r.) 


-14.8 ±1.0 ±24.8 


-10.3 ±0.6 ±19.6 




Ml (1-loop) 


-11.8 ± 1.1 ± 21.9 


-7.9 ±0.7 ± 17.3 




M2 (p.r.) 


-14.2 ±0.4 ±20.3 


-9.5 ± 0.3 ± 16.1 




M2 (1-loop) 


-13.6 ±0.6 ±24.3 


-8.5 ± 0.5 ± 18.0 




M3 (p.r.) 


-74.9 ±5.7 ±88.0 


-80.4 ±3.9 ±92.3 




M3 (1-loop) 


-72.6 ± 5.7 ± 142.0 


-73.5 ±4.0 ±131.0 





Table 5.4: e'/e results for Q 3 ensemble ((3 = 6.2). 



translates into a negative or very slightly positive value of e je (Tables |0] and pA 



and Fig. p.ll| ). 



Tables |5.3| and |5.4| contain our results for e'/e obtained in the above-described 



procedure. In all perturbative corrections we have used one-loop non-diagonal coef- 
ficients computed for gauge-noninvariant operators, which are assumed to be of the 
same order as those for gauge-invariant operators. The first quoted error is statistical 
(obtained by combining the individual errors in matrix elements by jackknife). The 
second error is the diagonal operator matching error due to uncertainty in the determi- 
nation of Zp and Zg. In order to estimate the non-diagonal operator matching error 
we compare two renormalization procedures: using strictly one-loop non-diagonal 
corrections (denoted "(1-loop)"), and resumming part of higher-order corrections in 
non-diagonal mixing by using non-perturbative renormalization factors Zp and Z$ 
(as explained in Section |5.3| ). The latter method is denoted "(p.r.)". 
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Finite volume and quenching effects were found small compared to noise. In 
addition to the quoted errors, there are uncertainties due to an unknown degree of 
validity of our ansatz in Sec. |5.3| , and due to unknown non-diagonal coefficients in the 
mixing matrix. The former error is uncontrolled at this point, since it is difficult to 
rigorously check our assumption made in Sec. The latter error is likely to be of 
the same order as the spread in e' / e between two different approaches to higher-order 



corrections (strictly one-loop and partial resummation), studied in Tables ^3| and fyA . 

Some other parameters used in obtaining these results are: ImAj = 1.5 • 1CT 4 , 
m t = 170 GeV, m b = 4.5 GeV, m c = 1.3 GeV, VL V+V , = 0.25, a^ =0) (2 GeV) = 0.195 
(the latter is based on setting the lattice scale by p meson mass). Short distance 
coefficients were obtained by two-loop running using the anomalous dimension and 



threshold matrices computed by Buras et al. [58 



Note that there are several choices to make in using our data to estimate e'/e. One 
can use the experimental values of KeA and ReA? in Eq. fl2.17| ) (Ml), or one can use 
the values obtained on the lattice (M3). One can also adopt an intermediate strategy 
(M2) of using the experimental amplitude ratio u and computed ReA^. When the 
higher-order chiral corrections are taken into account and the continuum limit is 
taken (so that u = 22), these three methods should converge. The method M2 is 
preferred, since in this method the overall error due to final state interactions cancels 
between real and imaginary parts of A amplitude, while the relative size of n and 
II2 contributions is given by the physical u>. 



Fig. |5.11| presents the data obtained on ensemble Qi, with error bars showing all 
three errors quoted in Tables and |5.4|, combined in quadrature. 
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Figure 5.11: A rough estimate of e'/e for Qi {(3 = 6.0) ensemble using the partially- 
nonperturbative procedure described in text. The Fermilab result's central value is 
shown with a horizontal line. The points are displaced horizontally for convenience. 
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Figure 5.12: A study of cutoff dependence of e'/e . Plotted are data obtained on 
(3 = 6.0 and (3 = 6.2 ensembles for M2 method. The error bars show only the 
statistical error in matrix elements and in Z$ and Zp constants. Systematic errors 
are significant but common to both ensembles. 



Fig. |5.12j contains comparison of results for e'/e on j3 = 6.0 and (3 = 6.2 ensembles. 
With current errors, the cutoff dependence of e'/e does not appear to be significant. 
The errors shown are purely statistical. 

5.5 Summary 

I have shown that perturbative operator matching conventionally employed in 
lattice calculations breaks down in the case of estimating e'/e . A partially nonper- 
turbative procedure is proposed as a temporary resolution. Light quark masses and 
(pseudo) scalar bilinear renormalization constants have been computed nonperturba- 
tively. 
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Our best result (at kaon mass, (3 = 6.0, M2 method) is 

Re(e'/e) = (-38.6 ± 9.3) • 10^ 4 . (5.25) 

The quoted error is a combination of the statistical and several known systematic 
errors. The reader should keep in mind that additional systematic errors present 
are very difficult to estimate. These errors are due to higher-order chiral terms (of 
the order of 40 %) and due to neglecting certain types of higher-order diagrams in 
partially nonperturbative operator matching. 

Taken at face value, our present estimate of e'/e would contradict the present 
experimental results, which favor a positive value. Therefore, the minimal Standard 
Model's description of direct CP violation would be proved to be invalid. At this 
point, it is too early to draw strong conclusions, since potentially large systematic 
uncertainties are present in our estimates. 

A recent study of e'/e with domain wall fermions by the RIKEN-BNL-Columbia 
group [|T^] has also produced negative values for e'/e . While this consistency is en- 



couraging (in particular, it shows that the assumptions made in our operator matching 
procedure do not change the qualitative picture), the reader should keep in mind that 
the errors due to using lowest-order chiral perturbation theory are common to both 
ours and RIKEN-BNL-Columbia calculations. 
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CHAPTER 6 



CONCLUSIONS 



In this thesis I report the techniques and results of work directed at calculating 
hadronic matrix elements of AS = 1 basis operators in Lattice QCD. Matrix elements 
of all basis operators have been computed with good statistical precision. Combined 
with lowest-order chiral perturbation theory, they allow one to make a theoretical 
prediction for AI = 1/2 rule and direct CP violation parameter e'/e . 

Based on the minimal Standard Model, we predict a substantial enhancement of 
AI = 1/2 transitions with respect to AI = 3/2 ones, consistent with experimental 
observations, although the exact amount of enhancement is subject to large system- 
atic uncertainties due to higher-order chiral terms. Thus our calculations confirm the 
currently dominant view that most of the enhancement comes from strong interac- 
tions. 

A theoretical prediction for e'/e is also made. Evaluated within the minimal Stan- 
dard Model framework, e'/e is found to have a negative value, in apparent conflict 
with experiments, although our prediction at present suffers from large systematic 
error. In the future, it is possible to improve the accuracy of this prediction by per- 
forming fully nonperturbative operator matching as well as including higher orders 
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in the chiral expansion. These improvements will come at a relatively easy compu- 
tational cost, since the main results of our calculation (matrix elements (7r + |Oj|A' + ) 
and (0\Oi\K°) in the lattice regularization scheme, which require most of the com- 
putational power) can be used along with above-mentioned improvements without 
repeating the calculations. If the sign of e'/e stays negative after these improve- 
ments, this would be the first failure of the minimal Standard Model to describe the 
real-world data, and so would have profound consequences in particle physics. 

In addition, we have computed light quark masses using nonperturbative renor- 
malization. 
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APPENDIX A 



EXPLICIT EXPRESSIONS FOR MATRIX ELEMENTS IN 
TERMS OF FERMION CONTRACTIONS. 



Operators in Eqs. (|2.2|) - |2.1ID can be decomposed into 1 = and 1 = 2 parts, 



which contribute, correspondingly, to A J = 1/2 and AI = 3/2 transitions. Here we 
give the expressions for these parts for completeness, since ReA , ReA 2 and e'/e are 

directly expressible in terms of their matrix elements. The 1 = parts are given as 
follows: 

Of = ^(s7 M (l-T5)rf)(^l-T5)w)-^(l-75)w)(w7 M (l-75)rf) 

+ I(^(l- 75 )d)(^(l- 75 )d) (A.l) 

Of = ^(^(l-75)n)(n 7 M (l-75)rf)-^(^(l-75)rf)(«7 M (l-75)n) 

+ l(s lfl (l- l5 )d)(d^(l- l5 )d) (A.2) 

of = (s % (l- l5 )d) E (g 7 ^(l-75)g) (A.3) 

q=u,d,s 

Of = (S a7/i (l-7*)d/j) E (^7 M (l-75)g a ) (A.4) 

q=u,d,s 

Of = (i 7M (l- 75 )d) 2 0?7 M (l+75)<?) (A.5) 

Of = (s a 7„(l-75)d/9) E (^7 M (l + 7 5 )g Q ) (A.6) 
Of = ^[(^(l-75)rf)(«7 M (l + 75)«)-(^(l- 75)^(^(1 + 75)^) 
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- (s7m(1-75M)(s7 m (1 + 75>)] (A.7) 
°8 0) = ^[0«7/.(l - 75)dp) (%7 M (1 + 7s)m q ) - (s a 7^(l - 7s)^) (%7 M (1 + 7s)d a ) 

- (s q7 ^(1- 75)^X^(1 + 75)^)] (A.8) 
Of = i[(Sr M (l - l 5 )d)(uY(l - 7b)u) - (St m (1 - 75)«)(U 7 ^(1 - 7s )d) 

- (^(1-75)^(^(1-75)5)] (A.9) 
Off = ^(^(l-75)«)(^(l-75)rf)-(^(l-75)rf)(«7 M (l-75)«) 

- (^(1-75)^(^(1-75)5)] (A.10) 

Expressions for 1 = 2 parts are as follows: 



of = 0? ) = |o 9 w = |oS = ^[(^(l-^)«)(^(l-7fe)d) 

+(s 7/ ,(l - 75)d)(w7 M (l - 75 )w) - (57a,(1 - 75)^)(rf7 M (l - 7s)<0] (A.ll) 

Of = i[(l7M(l-75)n)(U7 /1 (l + 75)rf) + (^(l-75)rf)(M7' 1 (l + 75)M) 

-(^(l-75)rf)(37 M (l + 75)rf)] (A.12) 

°8 = qlFaT/A 1 ~ 75)^/3)^/37^(1 + 75)^) + (S Q 7^(1 - 75)^/3) («/37 M (l + 7b)««) 

-(s a ^(l - 75)^(3^(1 + 75 )4)] (A.13) 
Of = Of = Of = Of = (A. 14) 

(Whenever the color indices are not shown, they are contracted within each bilinear, 
i.e. there are two color traces.) 



As mentioned in Sec. |3.2.3| , in order to compute matrix elements of I = op- 
erators one needs to evaluate three types of diagrams: "eight" (Fig. ^2|a), "eye" 
(Fig. |3.2| b) and "annihilation" (Fig. p.2| g). In the previous Appendix section we have 



given detailed expressions for computation of these contractions, given the spin-flavor 
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structure. Here we assign this structure to all contractions required for each operator, 
i.e. we express each matrix element in terms of contractions which were "built" in 
the previous section. 

Let us introduce some notation. The matrix element of the above operators have 
three components: 

2 2 

(tt+tt-IO^ ) = (Ei + Ii-S (2m«,) ) ±± (A.15) 

{Pn-pK)f 

where m is the common quark mass for s, d and u, and 

at = (A.16) 

Here E^ and Jj stand for "eight and "eye" contractions of the (n + \Oi\K + ) matrix 
element, Aj ~ (0\Oi\K°) / (m^ — m s ) is the "annihilation" diagram, S = (ir + \sd\K + ) 
is the "subtraction" diagram, and P = (0\s , -f5d\K ) is the two-point function. We 
compute a.i by averaging the ratio in the right-hand side of Eq. ([A. 16 ) over a suitable 
time range. 

Detailed expressions for Ei, Ii and are given below in terms of the basic con- 
tractions on the lattice. We label basic contractions by two letters, each representing 
a bilinear. For example, PP stands for contraction of the operator with spin structure 
( 75 ) (75), SS is for (1)(1), VV stands for ( 7/ J(7 M ), and AA is for (7^5) (7^5) • The 
staggered flavor is determined by the type of contraction, as explained in the previous 
Appendix section. Basic contractions are also labeled by their subscript. The first 
letter indicates whether it is an "eight" , "eye" or "annihilation" contraction, and the 
second is "U" for two, or "F" for one color trace. For example: PPeu stands for the 
"eight" contraction of the operator with spin-flavor structure (75 ® £5) (75 ® £5) with 
two color traces; VAaif stands for the "annihilation" contraction of the first type, in 
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which the derivative is taken with respect to quark mass on the external leg (see the 
previous Appendix section), the spin-flavor structure is (7^ <g> £5) (7^75 <8> 1), and one 
color trace is taken. What follows are the full expressions 6 . 
"Eight" parts: 



e[ 0) 


= \{VV EF + AA EF )-\{VV EU + AA EU ) 
3 3 


(A.17) 




= \{VV EU + AA EU )-\{W EF + AA EF ) 
3 3 


(A.18) 




= VVef + AA ef 


(A.19) 




= VV EU + AA EU 


(A.20) 


4 0) 


= 2{PP EF — SSef) 


(A.21) 




= 2{PP F u ~ SSeu) 


(A.22) 


Ei 0) 


= SSef - PPef + 1 -{VV EU - AA EU ) 


(A.23) 




= SSeu-PPeu + \{VVef-AA E f) 


(A.24) 




= -E§ } = ^(VVef + AAef -VV EU -AAeu) 


(A.25) 


E? 


= E? = ^E® = 2 -E { $ = 1 -{VV EU + AAeu + VV EF + AA EF ) 


(A.26) 


E { ? 


= E f ) = = E® = 


(A.27) 


E? 


= 1{AA EU -VV EU ) + SS EF -PP EF 


(A.28) 


Eg* 


= ^(AA EF -VV EF ) + SS EU -PPeu 


(A.29) 



"Eye" parts: 

/f) = Wnj + AAm (A.30) 
4 0) = VV IF + AA IF (A.31) 

6 Signs of operators O7 and 0$ have been changed in order to be consistent with the sign conven- 
tion of Buras et al. . 
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J< 0) 


= 3(VVrrr + AAtti) + 2(VVtp 




(A.32) 


/i 0) 


= 3(VVtt? + AAtw) + 2(Wttt 


+ AAtti) 


(A.33) 




= 2>(VVni - AAm) + AiPPrw 


- SStp) 


(A 34) 




= 3(Wrp - AArcO + AIPPttt 


- SSttt) 


(A.35) 




= UPPlT? - SStT?) 




(A.36) 




= 2(PP IU -SS IU ) 




(A.37) 


7 (o) 


= VV IF + AAj F 




(A.38) 


j(0) 
J 10 


= VV IU + AA IU 




(A.39) 



"Annihilation" parts are obtained by inserting the derivative with respect to 
— m s ) into every propagator involving the strange quark: 



4 0) 


= -{V A A1U + AV AIU ) 




(A.40) 


4 0) 


= -(VA a1f + AV A if) 




(A.41) 


4 0) 


= -3{VA AW + AVaw) ~ (VA A2U 


+ AV A2U ) 






-2(VA A1F + AVaif) ~ (VA A2F 


+ AV A2F ) 


(A.42) 


4 0) 


= -3{VA AlF + AV A1F ) - (VA A2F 


+ AV A2F ) 






-2(VA A1U + AV AW ) - (VA A2U 


+ AV A2U ) 


(A.43) 


4 0) 


= 3{VA AW - AVaw) + (VA A2U - 


AV A2U ) + 2(PS A2F 


- SP A2F )(A.U) 


4(0) 


= 3(VA A1F - AV A1F ) + (VA A2F - 


AV A2F ) + 2{PS A2U 


-SP A2U )(AA5) 


Af 


= \{VA A2U - AV A2U ) + {PS A2F - 


SP A2F ) 


(A.46) 


4 0) 


= \{VA A2F -AV A2F ) + {PS A2U - 


SP A2U ) 


(A.47) 


4 0) 


= VA A1F + AV A1F + ^ (VA A2U + AV A2U + VA A2F + AV A2F ) (A.48) 
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4? = VA A1U + AV A1U + ^(VA A2F + AV A2F + VA A2U + AV A2U ) (A.49) 
Of course, "eye" and "annihilation" contractions are not present in I = 2 operators. 
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